Underexpanded jet impingement in near vacuum environment

The investigation of the impact of a heated underexpanded exhaust plume impinging upon a flat surface in near-vacuum environments is of great significance owing to its relevance in the creation of solid landing pads on the Lunar surface within the next decade.¹ Retrorocket exhaust gases exert pressure and temperature loads on the surface, which can modify the landing pad and interact with the vehicle, causing structural and landing stability issues. The earliest experiments on highly underexpanded jet interactions with simulated Lunar surfaces were conducted by Stitt.² He studied cold air jets for a nozzle pressure ratio, i.e., ratio of stagnation pressure to background pressure, of up to 288 000. Vick et al.³ investigated the impact of highly underexpanded exhaust plumes at a Mach number of 5 on a perpendicular flat surface in the Langley 12.49 meter diameter vacuum sphere. He found that at a large stand-off distance, the surface shock appears to be a Mach disk and is closer to the nozzle exit than at a smaller stand-off distance, where the shock appears slightly above the flat surface. Additionally, he discovered that for low and medium nozzle stand-off distances, the surface pressure was directly proportional to the nozzle total pressure and unaffected by the ratio of nozzle total pressure to ambient pressure. At lower stand-off distances, the pressure distribution is a smooth bell curve. However, at larger stand-off distances, the distribution becomes oscillatory, with a pressure peak occurring at the outer edge of the plume interaction.³ 

Clark⁴ experimentally simulated the Apollo Lunar module rocket impingement on surfaces at very low stand-off heights. His results indicate the impingement pressure reaches a limiting maximum at the 2D_e stand-off distance, where D_e is the diameter of the nozzle exit, and the onset of nozzle choking at the 0.25D_e stand-off distance. The pressure distribution shape also changes from bell to annular below 2D_e. He also observed that maximum surface heating occurs at the 0.5D_e radial location rather than the nozzle axis. Inman et al.⁵ tested sonic and supersonic jet impingement at normal and oblique angles at various stand-off distances for nozzle pressure ratios of up to 634. The results indicate that in the presence of an impingement target, the onset location of flow unsteadiness can be diminished. It was also found that the flow instabilities grow more with a normal impingement than at oblique angles. Mehta et al.⁶ investigated underexpanded jet–surface interactions under Martian conditions. This analysis demonstrates that plumes released into the Martian environment provide significant surface pressure loads at high altitudes for the spacecraft, which is different from the plume–surface interaction which takes place under Lunar conditions. Vinze et al.⁷ found that the presence of the impingement shock significantly influences the local heat transfer over the plate for stand-off heights $≤$ 4D_e for nozzle pressure ratio ranges from 2.4 to 5.1. In an attempt to reduce the surface pressure on a flat surface in a rarefied environment, Kobayashi and Ukai⁸ used supersonic jet–jet interactions with an underexpanded nozzle flow and achieved a 48% reduction in stagnation pressure.

Alongside the experimental investigations, many researchers have conducted numerical analyses to examine the impingement effect under conditions of near vacuum. Roberts⁹ analytically investigated the dust dispersion phenomenon that resulted from plume impingement on the Lunar surface during the Apollo missions. Bird¹⁰ proposed a notable numerical simulation method called direct simulation Monte Carlo (DSMC) that can be used to solve non-continuum and non-equilibrium gas expansion in a near vacuum. Lumpkin et al.¹¹ used the DSMC technique to simulate the plume impingement of the Apollo Lunar module on the Lunar surface. Lane et al.¹² utilized both continuum computational fluid dynamics (CFD) and DSMC simulations of an Apollo-like Lunar lander to simulate the dynamics of regolith dust particle trajectories throughout the ascent and descent of lander rockets. Tosh et al.¹³ studied plume impingement on a flat surface to replicate the Lunar surface without taking erosion into account, using a hybrid continuum-rarefied CFD solver. To differentiate the rarefied flow from the continuum flow, Marichalar et al.¹⁴ employed a decoupled methodology utilizing DSMC and the Bird continuum breakdown parameter. Khasawneh et al.¹⁵ used DSMC simulations to verify the analysis of rarefied circular jet flow impinging on a flat plate. White et al.¹⁶ developed an uncoupled hybrid continuum-DSMC method and applied it to simulate the exhaust of a mono-propellant rocket positioned above the surface of an airless body.

Sharma et al.¹⁷ used a Navier–Stokes solver, together with the radiative transport equation (RTE), to analyze various phenomena such as engine–plume interactions and thermal loads on lander structures with clustered engine configurations, Mishra and Prasad¹⁸ studied disturbances on the Lunar surface induced by plume impingement, without taking the gas phase into account, using the Roberts erosion model. Gale et al.¹⁹ developed a simulation approach which integrates the impacts of shape irregularity of particles in a mixture using the discrete poly-disperse spherical model and the discrete element method (DEM). Rahimi et al.²⁰ examined the effect of the plume from a potential Lunar lander by applying physical conservation laws and using the Eulerian framework for the gas phase, as well as a particle-based discrete phase model (DPM) in the Lagrangian framework for the solid phase. Fontes et al.²¹ estimate the movement of Lunar regolith particles caused by contact with plumes using DSMC, considering various lander masses and altitudes. Ejtehadi et al.²² simulated the behavior of a two-phase flow using an Eulerian–Eulerian technique to anticipate how a rocket plume would interact with a dust blanket in a near-vacuum environment. Agir et al.²³ used the DSMC approach with the dsmcFoam+ solver to comprehensively investigate the impact of nozzle configuration and stagnation conditions on the interactions between plumes and surfaces. Cao et al. developed a solver called rarefiedMultiphaseFoam, which is based on dsmcFoam+,²⁴ implemented in OpenFOAM. They conducted tests on both steady and transient instances to evaluate their suitability for solving two-phase rarefied issues.²⁵ They then enhanced the solver to include solid–solid interactions and used it to analyze the interactions between a rocket exhaust plume and Lunar regolith.²⁶ 

Although numerical tools have greatly advanced in recent decades, experimental facilities remain in their early stages. One reason is that reproducing the extra-terrestrial conditions (reduced gravity, temperature, pressures, radiation environment, non-weathered regolith, etc.) in a terrestrial facility remains challenging.

Exhaust plumes from nozzles operating under near-vacuum conditions expand substantially and encounter a wide variety of flow regimes, from continuum flow inside the nozzle, the isentropic jet core adjacent to the nozzle exit, and the downstream regions that exhibit a change from continuum to free molecular flow. The different flow regimes of highly underexpanded flow into a vacuum and their flow impingement properties are shown in Fig. 1.

Due to the viscosity of the gas, a boundary layer develops along the nozzle walls and rapidly expands around the nozzle lip. Due to its high expansion, the boundary layer can cause a reverse flow behind the thruster, resulting in plume-induced contamination. The thickness of the boundary layer depends on the Reynolds number of the flow and is inversely proportional to the gas density. It is likely that in several instances, the backflow of the plume can come into contact with the surfaces of the spacecraft. An increased semi-divergence nozzle exit angle leads to a greater divergence of the plume boundary from the centerline. Ninety percent of the total mass flow is concentrated in a region of $30°$ from the nozzle axis for a conical nozzle with a 15° slant angle, but for a free jet from a sonic orifice, it is just 35%.²⁷ 

In a vacuum, the gas undergoes free expansion. This results in an increase in velocity and is accompanied by changes in temperature due to the redistribution of kinetic energy. The supersonic expansion downstream of the throat can create a compression shock system when the flow hits the nozzle wall. The shock goes into the freely expanding plume and changes the flow structure, especially the density and velocity fields. In reality, the flow is made up of both gas, unburned polymerized fuel droplets that were first sprayed into the combustion chamber, and droplets of combustion chamber cooling film, resulting in a two-phase non-uniform flow.

When a plume impinges on a flat surface in a vacuum, the gas transfers momentum and energy to the surface upon impact, leading to effects such as heating and potential erosion or modification of the surface material. The vacuum environment also influences heat transfer processes, affecting how generated heat is dissipated on the surface.

The Knudsen number, Kn, is one of the flow parameters used to determine the type of impingement. It is defined by the ratio between the mean free path of the gas to a reference length, i.e., $Kn=λL$⁠. The flow impingement is assumed to be free molecular when $Kn≥10$⁠, whereas it is assumed to be continuum when $Kn≤0.01$⁠. There is a transition flow impingement that occurs between these limits. When an object is placed within the continuum flow regime (⁠ $Kn≪1$⁠), a compression shock wave (surface shock) is generated in front of the object, see Fig. 1, and a layer of fluid is deflected along the surface of the object. A boundary layer develops on the surface and governs localized force and heat transfer. For the object situated in the free molecular regime, the gas molecules approaching the surface remain unaffected and undisturbed. Their contact with the wall involves the transfer of momentum and energy, as well as reflection or attachment, which is generally referred to as gas–surface interaction. When the value of Kn is reduced to less than 10, inter-molecular collisions become important, causing minor disruptions to the incoming flow.²⁷ 

The flow structure of rocket exhaust impingement quickly changes as the lander gets closer to the ground during a landing phase on an extraterrestrial body, such as the moon. When the lander is at a high altitude, the flow hitting the surface can be of low enough density that the flow is locally in the free molecular regime (⁠ $Kn≥10$⁠). As the lander descends toward the surface, the density of the flow impinging on the surface can become continuum (⁠ $Kn≤0.01$⁠). The present study aims to understand how the highly underexpanded exhaust plumes interacting with flat surfaces in near-vacuum environments affect surface pressure distributions and thermal loads at different heights.

The impingement study is timely, especially in light of the recent analytical study conducted by Metzger,^28,29 indicating that the volume of regolith liberated during the Apollo missions may have been far greater than previously estimated, emphasizing the current limitations in the knowledge and the need for conducting such a study to ensure more accurate forecasts of future landing operations.

Predicting underexpanded exhaust plume impingement is challenging due to limited experimental data in near-vacuum environments. The complexity arises from the restricted chamber capacity and the challenge of maintaining a consistently low back pressure during the whole plume discharge in the test facility. In this context, the present experimental work makes use of a large-volume, near-vacuum facility funded by the European Space Agency, specifically designed for conducting plume surface interaction (PSI) investigations in a near-vacuum environment. The facility, located within the Aerophysics Laboratory at the University of Glasgow, is comprised of many components, including a test chamber, a buffer chamber, connected conductance pipes, vacuum pumps, a control panel, and a dedicated software program. The schematic of the facility can be seen in Fig. 2. These components work together to regulate the essential valves and instruments required for the operation of the system and to maintain the vacuum level necessary for the experiment. The volume of the facility is separated into two main compartments: the test chamber, which is an axisymmetric chamber of 12 m³ volume with optical and electrical access ports, as shown in Fig. 3, and a cylindrical buffer tank combined with a conductance pipe, together creating a volume of 70 m³. Due to the restricted pumping capacity of the vacuum pumps, introducing mass during nozzle operation leads to an elevation in the background pressure inside the enclosed test chamber. Consequently, this has an impact on the structure of the plume within a short period and finally results in a scenario where the plume does not match the one seen in a vacuum environment. By including a buffer tank to increase the capacity, the test chamber's internal pressure is effectively controlled to ensure it remains within an acceptable range for a longer experiment duration, while minimizing any significant impact on the plume structure.^30,31 Additionally, as discussed in our previous work,^30,32 as long as the background pressure remains below 10 Pa as the mass is injected, the nozzle pressure ratio (NPR) remains greater than 100 000 and the flow can be assumed to be expanded in the same manner as it would under a pure vacuum. The above requirements resulted in the functioning of a supersonic nozzle for a duration of 2.5 seconds, in both constant and pulsed modes, with a mass flow rate of 1 g/s.

The current experiment involves the use of an underexpanded plume from a supersonic nozzle impinging on a flat surface in a near-vacuum environment. The schematic of the experimental setup is shown in Fig. 4. A total of eight experiments were performed in this study. The plume impingement is tested in steady and pulsed nozzle operation modes at stand-off distances of $y/De=4$ and $y/De=16$ and inclination angles of 0° and 30° for both stand-off distances. In steady mode, the nozzle is fired continuously for a duration of 2.5 seconds. In the pulsed mode, the solenoid operating frequency is configured to 6 Hz. The duration of each pulse is 0.083 seconds for both the flow-on and flow-off periods. The impingement plate is firmly fixed on an aluminum rig setup inside the chamber. The stand-off distance of the nozzle has been adjusted in two ways; for large movements, the whole support frame can be moved along the grooves of the aluminum rig support. For fine adjustments, a three-shaft traverse mechanism attached to the nozzle heater assembly can be adjusted. The angle of the plate is fixed using a high-precision electronic level (Bosch GLM 500 non-magnetic electronic meter). A laser alignment tool is used to align the nozzle center with the center of the plate. The final orientation of the nozzle and flat plate is ±0.3° in deviation from the requirement. The stand-off distance is measured from the exit of the nozzle to the center of the plate.

A scaled version of a simple conical supersonic convergent-divergent nozzle with fixed convergent-divergent angles and an exit Mach number of 6.6 was employed in this investigation. The diameter of the nozzle throat is 1 mm, while the exit diameter is 8.94 mm, with a tolerance of ±0.02 mm. Nitrogen gas is chosen to be the working gas for the nozzle due to its suitability in terms of safety. A pressurized nitrogen tank is connected to the nozzle entrance for this purpose through a stainless steel pipework assembly. Further in the line, a pressure regulator, a quick-open gate valve, and a mass flow meter have been attached. The nozzle discharge was initiated by a solenoid valve and a fast-acting double-actuated ball valve situated inside the vacuum chamber, which were both activated by a delay generator featuring a timed trigger.

The actuators are positioned in close proximity to the nozzle to mitigate any loss of total pressure caused by pipe wall friction and to minimize the pressure drop time when the valve is closed, which is particularly important for the pulsed mode operation. In addition to the pressure loss, the volume between the nozzle throat and ball valve significantly impacts the system's efficiency and response to changes in pressure and flow conditions, impacting the rise time—the time it takes from 10% (P10) to 90% (P90) of the system's fluid or gas to pass through. A larger volume results in longer P90 and P10 times due to increased buffer capacity, while a smaller volume leads to quicker response times. Flow dynamics are influenced by the pressure differential and flow rate through these components, with design considerations such as minimizing dead volumes and optimizing nozzle and valve sizes being essential for achieving desirable P90 and P10 times. A reduced pipe length and compact size of the heat exchanger resulted in a rapid increase in pressure with a calculated value of 11 milliseconds for the rise time. The rise time is within the specified range for pulsed operation and enough to activate pulse modulation mode.

To match the exit Reynolds number effectively under full-scale conditions, the temperature of the nitrogen gas is increased before entering the nozzle. To do this, a heat exchanger system is positioned between the gas supply pipe and the nozzle input. The heater system comprises a copper matrix shaft with parallel heat channels and a heater bundle. The matrix shaft is comprised of several heating channels, each with a diameter of 1 mm, to facilitate gas flow, as shown in Fig. 5. A heater bundle encloses the copper matrix to provide heat conduction. The heater bundle has cartridge heaters installed to evenly heat around the matrix shaft. A K-type thermocouple, linked to the heater bundle's outer surface monitors the temperature of the heater assembly. The feedback from the thermocouple was used to control the cartridge heaters through a proportional–integral–derivative (PID) controller (BriskHeat X2-120JTP model). The total temperature of the gas, T₀, is measured by another K-type thermocouple connected to the stagnation chamber of the nozzle as shown in Fig. 5. To minimize radiation losses, the heater bundle is insulated using a wool material. Using the heater, the nozzle is heated to a constant temperature of 900 K to obtain a Reynolds number more representative of a full-scale nozzle. The temperature of the heated gas is measured to be 750 K before each experimental run. This could be because heat is lost when it moves from the bundle to the heater matrix. It has been estimated³⁰ that the Reynolds number is greatly reduced from 38 668 for a gas temperature of 300 K to 9352 for a gas temperature of 900 K, hence, the temperature rise that occurs in this present study is thought to be sufficient to alter the flow Reynolds number.

The compact size of the heat exchanger will result in a rapid increase in pressure. Calculations based on the volume of the heater and nozzle assembly, lead to a rise time of 11 ms, which is within the specified range for pulsed operation (6 Hz) and enough to achieve pulse modulation mode. Furthermore, the gas temperature will be raised to 900 K at the heat exchanger, resulting in a pressure decrease of around 200 Pa. When the solenoid valve is activated, pressurized nitrogen gas flows into the heated copper matrix tubes of the heater assembly. The gas, at high temperature and pressure, proceeds into the stagnation chamber of the nozzle, resulting in a total pressure of 10 bar. The stagnation pressure is monitored through a Kulite pressure transducer attached to the opening of the stagnation chamber, as shown in Fig. 5.

A 45 cm diameter × 11 cm thick flat plate with provisions to install pressure transducers at different radial distances was used in this study. The material selected for the flat plate is PEEK (Poly Ether Ether Ketone), as its high emissivity allows for infrared images to be used to obtain temperature profiles on the flat plate. An image of the PEEK plate can be seen in Fig. 6. For pressure measurement, a Kulite pressure transducer is flush mounted at the center of the PEEK plate, which is aligned to the nozzle central axis, and another Kulite flush is mounted at 20 mm radially outward from the central point. Robert's analytical method³³ is used to choose the location of the outer transducer placement. Prior experience demonstrated that using the flat plate close to the hot nozzle could cause damage to the PEEK material due to radiative heat transfer. Hence, a millboard sheet (RS PRO FF700 Non-Ceramic Millboard of 3 mm thickness) cutout was placed beneath the heater, as depicted in Fig. 7, to impede the transfer of heat from the nozzle assembly to the PEEK plate.

Figure 8 shows the changes in chamber and nozzle properties during steady and pulsed mode operation. Prior to each experiment, the test chamber pressure is pumped down to near 1 Pa. Upon the release of mass during nozzle operation, the chamber pressure rises. A change in the gradient of the test chamber pressure increase was observed during the early stages of the plume release for both modes, as can be observed in Fig. 8(a). This is due to the difference between the pumped mass flow out of the chamber and the release of mass flow within the chamber. However, within a second, an equilibrium is reached, and the chamber pressure increases at a constant rate until it reaches close to 10 Pa for a steady 2.5 seconds of plume release and close to 8 Pa for 2 seconds of 6 Hz pulsing. Comparing the different modes, the steady operation introduces a greater amount of mass into the chamber compared with the pulsed operation, as evident in the temporal variations in pressure magnitude, as shown in the plots in Fig. 8(a). Note the pump speed is a function of pressure, but due to the tiny difference between modes in reaching the final chamber pressure (i.e., 8 and 10 Pa), we can safely assume that the pump operates at the same rate for both steady and pulsed modes. After the nozzle operation concludes, the chamber pressure continues to decrease as the pumps are active throughout.

A pressure regulator attached to the pipeline ensures 10 bar nozzle stagnation pressure for the entirety of the nozzle operation. In steady mode, the abrupt discharge of gas into the stagnation chamber due to the opening of the fast-acting ball valve caused the stagnation pressure to initially increase slightly over 10 bar, as seen in Fig. 8(b). The spike in pressure quickly subsides, and the pressure remains consistently at 10 bar until the nozzle shuts off at 2.5 seconds. During pulsed mode, the actuated fast-acting ball valve rapidly opens and closes, resulting in a variation in stagnation pressure ranging from 4 to 10 bar per cycle. The small overshoot note in steady mode can be detected in pulsed mode at the onset of every cycle. The gas temperature exhibited a marginal rise of about 10 K during the plume release, as seen in Fig. 8(c).

In real-world nozzle operation, the steady state is reached after a few seconds of nozzle firing. This happens because the nozzle wall is at the same temperature as the surrounding environment before combustion. Upon ignition, the nozzle wall temperature increases, affecting the transient behavior of the boundary layer and plume. When radiative heat loss from the walls and internal conduction of heat are balanced, a condition of thermal equilibrium is often reached. However, in the current experiment, the preheating of the nozzle allows the working gas to reach a stable condition soon after the firing [Fig. 8(c)].

The temporal variation of impingement pressure at 0° and 30° plate inclination at $y/De=4$ are presented in Figs. 9 and 10.

The impingement pressure measured on a flat plate at 0° inclination, i.e., the impingement plate is perpendicular to the plume axis, is shown in Fig. 9.

The pressure transducers are flush mounted at two radial distances; $x/De=0$ and $x/De=2.25$⁠. At $x/De=0$⁠, the transducers are exactly below the nozzle axis. A brief burst of high impingement pressure occurs concurrently with the ball valve opening, and the pressure quickly reaches nearly 5500 Pa for steady and pulsed modes, as seen in Fig. 9(a). It is typical for a starting jet to be accompanied by a moving primary shock wave that initiates the first contact, resulting in a sudden increase in pressure. An effort was made earlier in this work to perform Schlieren to detect the impingement of the shock wave. However, the extremely low background pressure and the decrease in pressure caused by the expansion of the plume render the shock wave invisible in the images. Nonetheless, the existence of a shock wave was detected by a hybrid CFD-DSMC simulation and will be elaborated upon in a subsequent section of this paper. Following the spike, there is a slow decline in impingement pressure noted in both modes during the initial second of nozzle operation, which then stabilizes for the remaining duration until the ball valve closes.

At $x/De=2.25$⁠, the flow tends to move radially outwards due to strong expansion fans from the highly underexpanded nozzle. The outward flow may move in a direction parallel to the surface of the second pressure sensor. In both modes, the second pressure sensor at $x/De=2.25$ measures a pressure of around 400 Pa during the initial moment of impingement. The spike observed at the centerline pressure is absent at the second sensor, as seen in Fig. 9(b). There is a consistent and gradual decrease in pressure as the measurement time advances. The direct cause of this occurrence is not apparent. However, it is hypothesized that this may be attributed to the interaction of high-temperature gas affecting the transducer in close proximity, resulting in a gradual decrease in pressure reading. Following the pulse, the pressure did not return to the background pressure but rather recorded a value below 0 Pa.

Figure 10 presents the results of pressure impingement measured at a plate inclined at 30° measured from the nozzle axis. The pressure sensors are again located at $x/De=0$ and $x/De=2.25$⁠. The initial transient peak is observed at the centerline pressure measurement in Fig. 10(a). Contrary to the pressure measurement at a 0° angle, the pressure does not decrease after the initial spike but quickly reaches a pseudo-steady state in both modes until the valve is completely closed. A marginal increase in pressure has been noted in the pulsed case. The magnitude decreased nearly 40% compared to the pressure measured at 0° inclination, see Fig. 10(a). A small increase in the magnitude of the pressure value is observed at $x/De=2.25$ in Fig. 10(b), however, the pressure profile is declining, as similarly observed for the pressure measurement in 0° inclination.

The impingement pressure results of 0° and 30° plate inclination from $y/De=16$ stand-off distance are presented in Figs. 11 and 12.

At higher stand-off distances, the centerline pressure readings for both pulsed and steady modes exhibit unique characteristics, as shown in Fig. 11(a). As the stand-off distance increases, the initial transient spike in pressure that was detected at a smaller stand-off distance is no longer present. This is likely because the impingement is occurring in the free molecular regime, where there is no existence of a shock wave to cause a pressure peak. The pressure sensors register a stable pressure of a magnitude near 500 Pa until the nozzle shuts off. It is interesting to note that in both modes, there was an abrupt peak in pressure near the end of the nozzle operation. The fact that the peak does not line up with the valve's activity is somewhat unexpected. Rather, for a short while, there was a pressure drop because of the lack of plume from the nozzle, but as the measurement came to a close, the pressure adversely increased. The cause of this pressure peak is uncertain.

At $x/De=2.25$ in Fig. 11(b), the pressure magnitude is slightly lower than the magnitude on 0° inclination at $y/De=4$⁠. While the pressure magnitude near the heated nozzle at $y/De=4$ consistently decreases, the same sensor does not exhibit the predicted heat impact as the stand-off distance is increased, as shown by the steady pressure for the entire duration of plume impingement. At the same location in the plot, the peak pressure that was detected at the end of the nozzle operation in centerline measurement may be seen in a very weak form.

The centerline measurement reveals a decrease in the magnitude of the surface pressure that has been seen for both modes in comparison to the flat plate at $y/De=4$⁠. The pressure peak that was seen at 0° inclination after the termination of the plume is also visible at 30° inclination. The measurement obtained at a 0° tilt has a somewhat higher magnitude than the value obtained at $x/De=2.25$⁠.

In this section, the results of an independent cold flow experiment conducted in steady firing mode at $y/De=4$ and 16 stand-off distances are compared in order to gain a clearer understanding of how heating the nozzle assembly affects the transducer readings. The experiment included deactivating the power supply to the heaters, which led to a gas temperature of 330 K. It is noteworthy to mention that the mass flow rate will increase, and the Reynolds number of the flow will be four times greater than that of the heated flow at this low temperature.³⁰ The transducers were positioned at $x/De=0$ radial distance as in the hot flow experiment.

The lack of heat causes the stagnation pressure to drop to a little less than 10 bar for the same supply pressure from the nitrogen tank. The impingement surface pressure value reflects the decrease in stagnation pressure, with the cold impingement pressure measuring somewhat lower than the hot impingement pressure. Figure 13(a) compares the centerline pressure between hot and cold flow. The cold flow exhibits an initial pressure increase that occurs when the ball valve opens, followed by steady pressure with a magnitude somewhat lower than the hot impingement results. A constant level of pressure is seen throughout the cold plume discharge, in contrast to the fluctuations seen in the hot flow. Despite the increase in stand-off to $y/De=16$⁠, the transducer reading may continue to exhibit the influence of the heated flow due to the marginal decrease in pressure over time, as illustrated in Fig. 13(b). The results of the cold flow impingement experiment imply that the temperature during the nozzle-fed hot gas discharge process affects the transducer reading. Therefore, to subsequently compare the experimental data with the steady-state analytical and computational results in the article, it was determined that an average of the impingement pressure readings from the beginning to the end of the plume discharge would be taken.

Infrared (IR) measurements on the impingement were recorded to understand the heat profile on the plate. Figure 14 compares the IR image of the flat plate before and after the plume impingement at the $y/De=16$ stand-off distance. The heated nozzle with heater wires can be seen at the top of the image, and the impingement plate is at the bottom.

The nozzle requires a few minutes to reach the desired temperature while using the heater. During this period of thermalization, the emission of radiation from the nozzle causes the flat plate surface temperature to increase. The distribution of heat before nozzle firing is visualized in Fig. 14(a) and the heat distribution during plume impingement is shown in Fig. 14(b).

The tilted flat plate positioned toward the IR camera enables a broader field of view in Fig. 15. In addition to the heat distribution, which is more pronounced at the upper portion of the plate due to its closer proximity, the infrared image also displays the reflection of the pressure ports that have been flushed with metal pins, as shown in Fig. 15(a). For tilted configuration at $y/De=16$⁠, the impingement of the hot plume creates an arc as shown in Fig. 15(b), possibly from the impingement of the outer ring of the tail shock.

The heat distribution taken from a line that goes from one end to the other through the center of the plate is compared in Figs. 16 and 17. This line is generated during post-processing in the FLIR software to compare the heat distribution along the plate surface.

The IR photos, Figs. 14 and 15, demonstrate that radiation from the heated nozzle assembly increases the flat plate temperature prior to plume impingement. The heat distribution before the plume release for a flat plate oriented at 0°, kept at $y/De=4$ stand-off height, resembles a bell curve that peaks at around 350 K and spreads to about $x/De=15$ radial distance on the surface, as shown in Fig. 16. The metal pills used to flush mount the transducer holes in the PEEK plate are the cause of the variation in background infrared radiation. With the plume impact, the limited area inside a $x/De=5$ radius receives additional heat, with a magnitude peaking at 400 K, as shown in Fig. 16. A decrease at the midpoint and a peak on both sides of the plot indicate the circular shape of the plume impingement. However, as the vertical distance increased to $y/De=16$⁠, a decrease in the maximum heat level was seen both before and after the plume impacts. In contrast to the bell curve distribution observed in $y/De=4$⁠, the higher stand-off gets a wider background heat distribution with a magnitude rising linearly from the plate edge and peaking at the center. The addition of plume heat increases the peak magnitude at the center point from 320 K to 330 K.

Due to the orientation of the plate in the tilted configuration, the top portion of the plate is positioned closer to the nozzle. This has resulted in an imbalance in the distribution of heat over the surface. Figure 17 plots the heat distribution along a horizontal line that passes through the middle of the circular plate at both stand-off distances. At $y/De=4$⁠, the profile of heat distribution and heat maximum for both background and plume impingement is similar to that of the profile in 0° inclination seen in Fig. 16, except for the fluctuation observed within $x/De=5$ radial distance. The fluctuations observed can be attributed to the existence of pressure slots within the specified region of the flat plate. At higher stand-off distances, the heat is distributed over the entire region of the circular plate, with the maximum heat registered at the center of the plate. The fluctuation in the background heat measurement is evident within $x/De=10$ radial distance for a higher stand-off distance, as shown in Fig. 17. The plate's orientation toward the IR camera places a strong emphasis on the slight emission variations brought about by the manufacturing roughness on the PEEK plate. Upon the impingement of a hot plume, the fluctuations in the measurement are more pronounced. A peak observed between $−10<x/De<0$ indicates the presence of impingement in an arc shape. The shape of the arc can be seen in Fig. 15.

An in-house, two-way coupled, hybrid CFD-DSMC solver³⁴ was used to simulate the $y/De=4$ stand-off height case at 0° nozzle inclination and with the steady firing mode. Comparisons of the obtained impingement pressures can then be made.

The hybrid simulation technique involves first simulating the entire domain using a CFD solver, in this case rhoCentralFoam. From this solution, the local Knudsen number, Kn, is used to determine the point at which continuum breakdown occurs within the flow field. The DSMC method, in this case dsmcFoam, is then used to simulate the rarefied region. For this specific case, in the final hybrid solution, CFD is used to calculate the flow field in the interior of the nozzle and the core of the jet, while the DSMC method is used for the rarefied regions beyond this.

First, an initial CFD simulation is conducted across the entirety of the domain. The domain is then decomposed following the aforementioned Kn condition into CFD and DSMC regions. A buffer region is also created, defined by the user through a chosen number of extension cells from the DSMC surface into the CFD region. Following this decomposition, the DSMC solver is applied across both the buffer and DSMC regions with the boundary conditions along the buffer-CFD interface set as fixed values of the macroscopic properties of temperature, pressure, and velocity that are obtained from the CFD solution. The entirety of the DSMC-buffer region is then initialized with DSMC particles from the fields of temperature, number density, velocity, heat flux, and stress from the initial CFD solution for the first iteration and the previous hybrid solution for any subsequent iterations. The DSMC simulation is then run until a steady state is reached and allowed to average over a sufficient number of time steps to reduce the statistical scatter in the final macroscopic fields. Being two-way coupled, this DSMC solution is then used to initialize the CFD solver across the entirety of the CFD-buffer region in the same manner. The boundary conditions for the DSMC-buffer interface are set as fixed values of the macroscopic properties of temperature, pressure, and velocity that are obtained from the DSMC simulation. The CFD simulation is then run until a steady state is obtained. This algorithm of domain decomposition, DSMC, and CFD is then repeated until the CFD and DSMC solutions converge within the buffer region and a final steady-state hybrid solution is thus obtained. For a more rigorous description of the hybrid CFD-DSMC solver, the reader can consult the work of Vasileiadis and White.³⁴ 

An axisymmetric simulation was performed around the nozzle centerline. The continuum breakdown was selected to be where Kn reached 0.05 or greater, $Kn≥0.05$⁠, and the characteristic length, L, to calculate the mean free path-based Kn was the nozzle throat diameter of 1 mm. A buffer layer ten cells wide was used. The boundary conditions for the initial CFD simulation are set in accordance with the conducted experiments; the stagnation temperature and pressure at the inlet of the nozzle are set to 780 K and 10 bar, respectively. The nozzle walls are prescribed a temperature of 900 K, a Maxwell slip condition for velocity, and a Smoluchowski jump condition for temperature. The impingement plate is also prescribed a Maxwell slip condition for velocity and a zero gradient condition for both pressure and temperature. The ambient boundary is assigned a wave transmissive condition for pressure and temperature and an inlet-outlet condition for velocity. A schematic of the domain and the assigned boundary types can be seen in Fig. 18.

For the DSMC simulation, the nozzle walls and impingement plate are set as a diffuse wall at the same temperatures as in the CFD. Each DSMC particle represented 1.25  $×1010$ real nitrogen molecules and cell-weighting ensured that each cell initially contained 20 DSMC particles. A time step of $1×10−9$ s was used. Two rotational and two vibrational degrees of freedom for nitrogen were used in both the DSMC and CFD simulations. A mean free path-based automatic mesh refinement was used, and the final refined hybrid mesh was comprised of 115 247 cells (550 prismatic elements along the axis of symmetry, 2177 polyhedral elements in regions of refinement, and 112 520 hexahedral elements). This final hybrid mesh can be seen in Fig. 20.

Contours of velocity magnitude, streamlines, and Mach number are presented in Figs. 21, 22, and 23, respectively. The exit Mach number in the simulation case is 6.1, lower than the theoretical design Mach number of 6.6. This is possibly due to the fact that the boundary layer encompasses a relatively large proportion of the nozzle's small cross-sectional area (most clearly seen in the Mach contour in Fig. 23), resulting in a reduction in nozzle performance.³⁶ There is no measurement of the exit Mach number in the experiment so a comparison cannot be made to investigate this possibility.

The underexpanded plume expands rapidly beyond the nozzle exit and reaches approximately Mach 10 upstream of the stand-off or bow shock, the presence of which was implied by the experimental measurements.

As mentioned previously, the locations where the flow becomes rarefied are determined by the local Knudsen number exceeding a value of 0.05, $Kn≥0.05$⁠. A diagram of the domain showing these calculated regions is shown in Fig. 24. For almost the entirety of the interior region of the nozzle, and the core of the jet down to the point it impinges on the flat plate, the flow is in the continuum region.^37,38 The flow becomes rarefied as it is expanding outwards radially. There is also a small degree of rarefaction along the interior of the nozzle wall toward the nozzle exit.

A comparison of the impingement pressure profiles on the surface of the plate obtained through simulation, experiment, and the aforementioned analytical method of Roberts is made in Fig. 25. All three of the methods are in relatively good agreement. Measurements of 4849 and 4900 Pa were obtained for the centerline pressure in the simulation and experiment, respectively, which is a difference of 0.7%. At the location of the second experimental probe, a pressure of 370 Pa was measured, whilst the simulation predicted a pressure of 250 Pa. These measurements are compared with the Roberts method with exit Mach numbers of both 6.6 and 6.1, corresponding to the theoretical exit Mach number and that which was obtained in the simulation case. It is worth noting that the Roberts method assumes the expansion of the plume is into a perfect vacuum, which differs from both the experiment and simulation cases where a background pressure of 1.8 Pa, which increases throughout the experiment, is present.

Due to the optical thickness of the gas being very low under the tested lunar conditions, it was not possible to conduct Schlieren measurements for the plume–surface interaction experimentally. As the obtained impingement pressures of the experiment and simulation were in good agreement, the numerical results can be used to generate a numerical Schlieren image to provide an indication of the shock structure present during the interaction. The case of $0°$ inclination and $4De$ stand-off height was studied in this manner.

The formation of a bow shock is once again confirmed by the numerical Schlieren in Fig. 26. The shock is formed at a height of approximately 2 mm above the surface of the flat plate and is approximately 100 mm in width. Furthermore, it can be seen that a barrel shock is also present within the plume structure. The maximum width of the barrel shock is approximately 8 mm and occurs approximately 2.5 mm downstream of the nozzle exit plane.

When a rocket exhaust expands into a rarefied environment, the pressure profile and structure of a PSI change from those in a terrestrial environment. Moreover, the direct impingement of a thruster plume on the surface of a spacecraft is responsible for the emergence of undesired torques, localized surface heating, and surface contamination. Hence, understanding the PSI effects under near-vacuum environmental conditions is crucial for future lander missions to celestial bodies with no or thin atmospheres. The present study investigates the complex flow patterns resulting from underexpanded axisymmetric plumes that impinge on flat plates in a near vacuum environment. The experiments were conducted within a specialized plume regolith facility at the University of Glasgow. A flat surface made up of PEEK equipped with pressure transducers is installed at different stand-off heights (4 and 16  $y/De$⁠) and inclination (0° and 30°) below a scaled version of a convergent-divergent nozzle. The nozzle stagnation chamber is heated to attain Reynolds number similitude with a full-scale nozzle. The temperature profile of the plume impingement on the plate surface was recorded using an infrared camera.

An initial chamber pressure of near 1 Pa increased to 10 Pa at the end of 2.5 seconds of steady nozzle firing at a mass flow rate of 1 g/s but with intermittent mass addition in pulsed firing, the chamber pressure saw an increase up to 7 Pa. A high NPR of 100 000 ensures that the flow resembles that in the pure vacuum at a specified stand-off height. A quasi-steady impingement pressure up to 5500 Pa was registered for both steady and pulsed modes at lower stand-off height. The presence of a primary shock wave is indicated by an initial increase in the pressure at the centerline of the lower stand-off measurement. This shock wave is identified in the numerical rebuild results, but cannot be seen in the Schlieren visualization owing to the low density of the surrounding environment. With the flow moving radially outward, the pressure at $x/De=2.25$ is 10 times lower than the centerline pressure. The quasi-steadiness in the temporal pressure measurement may be attributed to the interaction of a hot plume with the pressure sensor. A separate cold flow experiment has identified the impact of high-temperature gas on pressure measurement. The centerline pressure decreased nearly to 40% when the flat plate was inclined to 30°. At the higher stand-off distance, the initial spike due to the primary shock wave is absent, but an abrupt peak at the end of nozzle firing has been registered at centerline measurement in both plate inclinations.

The emission of infrared radiation from the flat plate indicates the existence of heat prior to the firing of the nozzle, which is caused by the release of heat radiation from the heated nozzle. A brighter center region indicates the impingement of a hot plume on the PEEK plate. When the plate is inclined, the interplay of shock waves may cause the plume to generate an arc on the top section of the plate. The temperature distribution along the diameter of the circular plate follows a bell-shaped curve, with the highest value in the center before the hot plume hits it. However, after the hot plume impinges on the plate, there is a decrease in temperature at the midway and peaks on both sides, indicating the presence of a shock wave at a lower stand-off distance. Compared to the bell curve distribution in $y/De=4$⁠, the greater stand-off results in a larger background heat distribution with a linear increase in amplitude.

An in-house hybrid CFD-DSMC solver was used to model the $y/De=4$ stand-off height scenario at 0° nozzle inclination and steady firing mode, validating the experimental results. The hybrid simulation methodology employs the Knudsen number, Kn, to identify continuum breakdown points in the flow field. It includes a Navier–Stokes–Fourier equation solver modeling of the continuum zone and the DSMC method simulating the rarefied region. The computational analysis indicates that the Mach number exceeds 10 beyond the nozzle exit, and a surface shock forms above the flat plate, as indicated by the experimental data. The impingement pressure profile from computational and experimental results is then compared with an analytical method proposed by Roberts, which assumes plume growth into a perfect vacuum, unlike the experiment and simulation scenarios, which include a background pressure of 1.8 Pa that rises during the experiment. It was shown that all three methods are in relatively good agreement.

The authors wish to express gratitude toward the European Space Agency (ESA) for funding the experimental part of this research through ESA Contract No. 4000115469/15/NL/KML/fg, Engineering and Physical Sciences Research Council for funding the computational work through EPSRC Grant number EP/V012010/1, BUSCH UK for their ongoing assistance in ensuring the efficient functioning of the facility, Dr. Muhammed Burak Agir and Dr. H. Zare-Behtash, for their technical and operational support in the initial stages of the experimental work, Ph.D. students of aero-physics laboratory (Ting Tsung Chang and Gaargi Jain), and the technicians of the University of Glasgow James Watt School of Engineering.

The authors have no conflicts to disclose.

S. Subramanian: Investigation (lead); Writing – original draft (lead). B. Craig: Investigation (supporting); Validation (lead); Writing – original draft (equal); Writing – review & editing (supporting). C. White: Investigation (supporting); Writing – review & editing (lead). K. Kontis: Supervision (lead); Writing – review & editing (supporting). D. Evans: Investigation (supporting). J. Van den Eynde: Investigation (supporting); Project administration (lead); Writing – review & editing (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.