This dissertation develops and demonstrates a new physics-based approach that provides computational stability of subgrid stress models in large eddy simulations while producing far smaller changes in the original subgrid stress and subgrid production fields than do current \textit{ad hoc} stabilization methods. A pseudo-spectral code that is shown here to be almost entirely non-dissipative yet inherently stable without any subgrid model is used to conduct simulations with stable and unstable subgrid stress models. Results show that initial instability, subsequent exponential growth, and eventual machine overflow occur via a highly localized dynamical process that results from interactions among terms in the kinetic energy and enstrophy transport equations. This process begins first at one material point and then occurs at increasingly more material points, with local exponential growth rates of kinetic energy and enstrophy being the same for all points, until machine overflow eventually occurs at the material point where the process began first. A Lagrangian backtracking scheme is developed and applied to this material point, allowing backward-in-time tracking of all terms in the kinetic energy and enstrophy transport equations. This gives insights into the dynamics that produce this local instability and its subsequent exponential growth, with the initial instability shown to result from interactions between the subgrid production and subgrid redistribution terms. Elementary backscatter limiting based on locally reducing individual subgrid stress components that contribute to local kinetic energy backscatter is shown to stabilize any stress model, but still produces substantial changes in the stress and production fields. The rational Boolean stabilization method developed here instead uses the local subgrid production and subgrid redistribution rates to determine where and how individual subgrid stress components must be rescaled to provide local backscatter limiting and/or forward scatter amplification. This stabilizes all subgrid stress models while producing only small changes in the subgrid stress and production fields. Rational Boolean stabilization is computationally fast, and can be generalized to stabilize models for other subgrid terms in large eddy simulations while producing only small changes in their resulting fields. This solves a key problem that has previously limited the accuracy of large eddy simulations.

The Magnetoplasmadynamic (MPD) thruster is an electromagnetic thruster that produces a higher specific impulse than conventional chemical rockets and greater thrust densities than electrostatic thrusters, but the well-known operational limit---referred to as ``onset"---imposes a severe limitation efficiency and lifetime. This phenomenon is associated with large fluctuations in operating voltage, high rates of electrode erosion, and three-dimensional instabilities in the plasma flow-field which cannot be adequately represented by two-dimensional, axisymmetric models. Simulations of the Princeton Benchmark Thruster (PBT) were conducted using the three-dimensional version of the magnetohydrodynamic (MHD) code, MACH. Validation of the numerical model is partially achieved by comparison to equivalent simulations conducted using the well-established two-dimensional, axisymmetric version of MACH. Comparisons with available experimental data was subsequently performed to further validate the model and gain insights into the physical processes of MPD acceleration. Thrust, plasma voltage, and plasma flow-field predictions were calculated for the PBT operating with applied currents in the range $6.5kA < J < 23.25kA$ and mass-flow rates of $1g/s$, $3g/s$, and $6g/s$. Comparisons of performance characteristics between the two versions of the code show excellent agreement, indicating that MACH3 can be expected to be as predictive as MACH2 has demonstrated over multiple applications to MPD thrusters. Predicted thrust for operating conditions within the range which exhibited no symptoms of the onset phenomenon experimentally also showed agreement between MACH3 and experiment well within the experimental uncertainty. At operating conditions beyond such values , however, there is a discrepancy---up to $\sim20\%$---which implies that certain significant physical processes associated with onset are not currently being modeled. Such processes are also evident in the experimental total voltage data, as is evident by the characteristic ``voltage hash", but not present in predicted plasma voltage. Additionally, analysis of the predicted plasma flow-field shows no breakdown in azimuthal symmetry, which is expected to be associated with onset. This implies that perhaps certain physical processes are modeled by neither MACH2 nor MACH3; the latter indicating that such phenomenon may not be inherently three dimensional and related to the plasma---as suggested by other efforts---but rather a consequence of electrode material processes which have not been incorporated into the current models.