General-domain compressible Navier-Stokes solvers exhibiting quasi-unconditional stability and high-order accuracy in space and time

This thesis presents a new class of solvers for the subsonic compressible Navier-Stokes equations in general two- and three-dimensional spatial domains. The proposed methodology incorporates: 1) A novel linear-cost implicit solver based on use of higher-order backward differentiation formulae (BDF) and the alternating direction implicit approach (ADI); 2) A fast explicit solver; 3) Dispersionless spectral spatial discretizations; and 4) A domain decomposition strategy that negotiates the interactions between the implicit and explicit domains. In particular, the implicit methodology is quasi-unconditionally stable (it does not suffer from CFL constraints for adequately resolved flows), and it can deliver orders of time accuracy between two and six in the presence of general boundary conditions. In fact this thesis presents, for the first time in the literature, high-order time-convergence curves for Navier-Stokes solvers based on the ADI strategy---previous ADI solvers for the Navier-Stokes equations have not demonstrated orders of temporal accuracy higher than one. An extended discussion is presented in this thesis which places on a solid theoretical basis the observed quasi-unconditional stability of the methods of orders two through six. The performance of the proposed solvers is favorable. For example, a two-dimensional rough-surface configuration including boundary layer effects at Reynolds number equal to one million and Mach number 0.85 (with a well-resolved boundary layer, run up to a sufficiently long time that single vortices travel the entire spatial extent of the domain, and with spatial mesh sizes near the wall of the order of one hundred-thousandth the length of the domain) was successfully tackled in a relatively short (approximately thirty-hour) single-core run; for such discretizations an explicit solver would require truly prohibitive computing times. As demonstrated via a variety of numerical experiments in two- and three-dimensions, further, the proposed multi-domain parallel implicit-explicit implementations exhibit high-order convergence in space and time, useful stability properties, limited dispersion, and high parallel efficiency.

I. Stokes flow past a thin screen. II. Viscous flows past porous bodies of finite size

Part I

The slow, viscous flow past a thin screen is analyzed based on Stokes equations. The problem is reduced to an associated electric potential problem as introduced by Roscoe. Alternatively, the problem is formulated in terms of a Stokeslet distribution, which turns out to be equivalent to the first approach.

Special interest is directed towards the solution of the Stokes flow past a circular annulus. A "Stokeslet" formulation is used in this analysis. The problem is finally reduced to solving a Fredholm integral equation of the second kind. Numerical data for the drag coefficient and the mean velocity through the hole of the annulus are obtained.

Stokes flow past a circular screen with numerous holes is also attempted by assuming a set of approximate boundary conditions. An "electric potential" formulation is used, and the problem is also reduced to solving a Fredholm integral equation of the second kind. Drag coefficient and mean velocity through the screen are computed.

Part II

The purpose of this investigation is to formulate correctly a set of boundary conditions to be prescribed at the interface between a viscous flow region and a porous medium so that the problem of a viscous flow past a porous body can be solved.

General macroscopic equations of motion for flow through porous media are first derived by averaging Stokes equations over a volume element of the medium. These equations, including viscous stresses for the description, are more general than Darcy's law. They reduce to Darcy's law when the Darcy number becomes extremely small.

The interface boundary conditions of the first kind are then formulated with respect to the general macroscopic equations applied within the porous region. An application of such equations and boundary conditions to a Poiseuille shear flow problem demonstrates that there usually exists a thin interface layer immediately inside the porous medium in which the tangential velocity varies exponentially and Darcy's law does not apply.

With Darcy's law assumed within the porous region, interface boundary conditions of the second kind are established which relate the flow variables across the interface layer. The primary feature is a jump condition on the tangential velocity, which is found to be directly proportional to the normal gradient of the tangential velocity immediately outside the porous medium. This is in agreement with the experimental results of Beavers, et al.

The derived boundary conditions are applied in the solutions of two other problems: (1) Viscous flow between a rotating solid cylinder and a stationary porous cylinder, and (2) Stokes flow past a porous sphere.

The flow between rotating coaxial disks

Numerical approximations of nonunique solutions of the Navier-Stokes equations are obtained for steady viscous incompressible axisymmetric flow between two infinite rotating coaxial disks. For example, nineteen solutions have been found for the case when the disks are rotating with the same speed but in opposite direction. Bifurcation and perturbed bifurcation phenomena are observed. An efficient method is used to compute solution branches. The stability of solutions is analyzed. The rate of convergence of Newton's method at singular points is discussed. In particular, recovery of quadratic convergence at "normal limit points" and bifurcation points is indicated. Analytical construction of some of the computed solutions using singular perturbation techniques is discussed.

Some problems of edge waves and standing waves on beaches

Some problems of edge waves and standing waves on beaches are examined.

The nonlinear interaction of a wave normally incident on a sloping beach with a subharmonic edge wave is studied. A two-timing expansion is used in the full nonlinear theory to obtain the modulation equations which describe the evolution of the waves. It is shown how large amplitude edge waves are produced; and the results of the theory are compared with some recent laboratory experiments.

Traveling edge waves are considered in two situations. First, the full linear theory is examined to find the finite depth effect on the edge waves produced by a moving pressure disturbance. In the second situation, a Stokes' expansion is used to discuss the nonlinear effects in shallow water edge waves traveling over a bottom of arbitrary shape. The results are compared with the ones of the full theory for a uniformly sloping bottom.

The finite amplitude effects for waves incident on a sloping beach, with perfect reflection, are considered. A Stokes' expansion is used in the full nonlinear theory to find the corrections to the dispersion relation for the cases of normal and oblique incidence.

Finally, an abstract formulation of the linear water waves problem is given in terms of a self adjoint but nonlocal operator. The appropriate spectral representations are developed for two particular cases.

Compressible flows at small Reynolds numbers

The problem of the slow viscous flow of a gas past a sphere is considered. The fluid cannot be treated incompressible in the limit when the Reynolds number Re, and the Mach number M, tend to zero in such a way that Re ~ o(M^2 ). In this case, the lowest order approximation to the steady Navier-Stokes equations of motion leads to a paradox discovered by Lagerstrom and Chester. This paradox is resolved within the framework of continuum mechanics using the classical slip condition and an iteration scheme that takes into account certain terms in the full Navier-Stokes equations that drop out in the approximation used by the above authors. It is found however that the drag predicted by the theory does not agree with R. A. Millikan's classic experiments on sphere drag.

The whole question of the applicability of the Navier-Stokes theory when the Knudsen number M/Re is not small is examined. A new slip condition is proposed. The idea that the Navier-Stokes equations coupled with this condition may adequately describe small Reynolds number flows when the Knudsen number is not too large is looked at in some detail. First, a general discussion of asymptotic solutions of the equations for all such flows is given. The theory is then applied to several concrete problems of fluid motion. The deductions from this theory appear to interpret and summarize the results of Millikan over a much wider range of Knudsen numbers (almost up to the free molecular or kinetic limit) than hitherto Believed possible by a purely continuum theory. Further experimental tests are suggested and certain interesting applications to the theory of dilute suspensions in gases are noted. Some of the questions raised in the main body of the work are explored further in the appendices.

Nonexistence of looping trajectories in hydromagnetic waves of finite amplitude. Breaking of waves in a cold collision-free plasma in a magnetic field. On stabilities of periodic waves in a cold collision-free plasma in a magnetic field.

This dissertation consists of three parts. In Part I, it is shown that looping trajectories cannot exist in finite amplitude stationary hydromagnetic waves propagating across a magnetic field in a quasi-neutral cold collision-free plasma. In Part II, time-dependent solutions in series expansion are presented for the magnetic piston problem, which describes waves propagating into a quasi-neutral cold collision-free plasma, ensuing from magnetic disturbances on the boundary of the plasma. The expansion is equivalent to Picard's successive approximations. It is then shown that orbit crossings of plasma particles occur on the boundary for strong disturbances and inside the plasma for weak disturbances. In Part III, the existence of periodic waves propagating at an arbitrary angle to the magnetic field in a plasma is demonstrated by Stokes expansions in amplitude. Then stability analysis is made for such periodic waves with respect to side-band frequency disturbances. It is shown that waves of slow mode are unstable whereas waves of fast mode are stable if the frequency is below the cutoff frequency. The cutoff frequency depends on the propagation angle. For longitudinal propagation the cutoff frequency is equal to one-fourth of the electron's gyrofrequency. For transverse propagation the cutoff frequency is so high that waves of all frequencies are stable.

Theory of viscous and thermal attenuation of sound by small spheres

The problem is to calculate the attenuation of plane sound waves passing through a viscous, heat-conducting fluid containing small spherical inhomogeneities. The attenuation is calculated by evaluating the rate of increase of entropy caused by two irreversible processes: (1) the mechanical work done by the viscous stresses in the presence of velocity gradients, and (2) the flow of heat down the thermal gradients. The method is first applied to a homogeneous fluid with no spheres and shown to give the classical Stokes-Kirchhoff expressions. The method is then used to calculate the additional viscous and thermal attenuation when small spheres are present. The viscous attenuation agrees with Epstein's result obtained in 1941 for a non-heat-conducting fluid. The thermal attenuation is found to be similar in form to the viscous attenuation and, for gases, of comparable magnitude. The general results are applied to the case of water drops in air and air bubbles in water.

For water drops in air the viscous and thermal attenuations are camparable; the thermal losses occur almost entirely in the air, the thermal dissipation in the water being negligible. The theoretical values are compared with Knudsen's experimental data for fogs and found to agree in order of magnitude and dependence on frequency. For air bubbles in water the viscous losses are negligible and the calculated attenuation is almost completely due to thermal losses occurring in the air inside the bubbles, the thermal dissipation in the water being relatively small. (These results apply only to non-resonant bubbles whose radius changes but slightly during the acoustic cycle.)

The Equation of State and Petrogenesis of Komatiite

(1) Equation of State of Komatiite

The equation of state (EOS) of a molten komatiite (27 wt% MgO) was detennined in the 5 to 36 GPa pressure range via shock wave compression from 1550°C and 0 bar. Shock wave velocity, U_S, and particle velocity, U_P, in km/s follow the linear relationship U_S = 3.13(±0.03) + 1.47(±0.03) U_P. Based on a calculated density at 1550°C, 0 bar of 2.745±0.005 glee, this U_S-U_P relationship gives the isentropic bulk modulus K_S = 27.0 ± 0.6 GPa, and its first and second isentropic pressure derivatives, K'_S = 4.9 ± 0.1 and K"_S = -0.109 ± 0.003 GPa^-1.

The calculated liquidus compression curve agrees within error with the static compression results of Agee and Walker [1988a] to 6 GPa. We detennine that olivine (FO₉₄) will be neutrally buoyant in komatiitic melt of the composition we studied near 8.2 GPa. Clinopyroxene would also be neutrally buoyant near this pressure. Liquidus garnet-majorite may be less dense than this komatiitic liquid in the 20-24 GPa interval, however pyropic-garnet and perovskite phases are denser than this komatiitic liquid in their respective liquidus pressure intervals to 36 GPa. Liquidus perovskite may be neutrally buoyant near 70 GPa.

At 40 GPa, the density of shock-compressed molten komatiite would be approximately equal to the calculated density of an equivalent mixture of dense solid oxide components. This observation supports the model of Rigden et al. [1989] for compressibilities of liquid oxide components. Using their theoretical EOS for liquid forsterite and fayalite, we calculate the densities of a spectrum of melts from basaltic through peridotitic that are related to the experimentally studied komatiitic liquid by addition or subtraction of olivine. At low pressure, olivine fractionation lowers the density of basic magmas, but above 14 GPa this trend is reversed. All of these basic to ultrabasic liquids are predicted to have similar densities at 14 GPa, and this density is approximately equal to the bulk (PREM) mantle. This suggests that melts derived from a peridotitic mantle may be inhibited from ascending from depths greater than 400 km.

The EOS of ultrabasic magmas was used to model adiabatic melting in a peridotitic mantle. If komatiites are formed by >15% partial melting of a peridotitic mantle, then komatiites generated by adiabatic melting come from source regions in the lower transition zone (≈500-670 km) or the lower mantle (>670 km). The great depth of incipient melting implied by this model, and the melt density constraint mentioned above, suggest that komatiitic volcanism may be gravitationally hindered. Although komatiitic magmas are thought to separate from their coexisting crystals at a temperature =200°C greater than that for modern MORBs, their ultimate sources are predicted to be diapirs that, if adiabatically decompressed from initially solid mantle, were more than 700°C hotter than the sources of MORBs and derived from great depth.

We considered the evolution of an initially molten mantle, i.e., a magma ocean. Our model considers the thermal structure of the magma ocean, density constraints on crystal segregation, and approximate phase relationships for a nominally chondritic mantle. Crystallization will begin at the core-mantle boundary. Perovskite buoyancy at > 70 GPa may lead to a compositionally stratified lower mantle with iron-enriched mangesiowiistite content increasing with depth. The upper mantle may be depleted in perovskite components. Olivine neutral buoyancy may lead to the formation of a dunite septum in the upper mantle, partitioning the ocean into upper and lower reservoirs, but this septum must be permeable.

(2) Viscosity Measurement with Shock Waves

We have examined in detail the analytical method for measuring shear viscosity from the decay of perturbations on a corrugated shock front The relevance of initial conditions, finite shock amplitude, bulk viscosity, and the sensitivity of the measurements to the shock boundary conditions are discussed. The validity of the viscous perturbation approach is examined by numerically solving the second-order Navier-Stokes equations. These numerical experiments indicate that shock instabilities may occur even when the Kontorovich-D'yakov stability criteria are satisfied. The experimental results for water at 15 GPa are discussed, and it is suggested that the large effective viscosity determined by this method may reflect the existence of ice VII on the Rayleigh path of the Hugoniot This interpretation reconciles the experimental results with estimates and measurements obtained by other means, and is consistent with the relationship of the Hugoniot with the phase diagram for water. Sound waves are generated at 4.8 MHz at in the water experiments at 15 GPa. The existence of anelastic absorption modes near this frequency would also lead to large effective viscosity estimates.

(3) Equation of State of Molybdenum at 1400°C

Shock compression data to 96 GPa for pure molybdenum, initially heated to 1400°C, are presented. Finite strain analysis of the data gives a bulk modulus at 1400°C, K'_S. of 244±2 GPa and its pressure derivative, K'_OS of 4. A fit of shock velocity to particle velocity gives the coefficients of U_S = C_O+S U_P to be C_O = 4.77±0.06 km/s and S = 1.43±0.05. From the zero pressure sound speed, C_O, a bulk modulus of 232±6 GPa is calculated that is consistent with extrapolation of ultrasonic elasticity measurements. The temperature derivative of the bulk modulus at zero pressure, θK_OSθT|_P, is approximately -0.012 GPa/K. A thermodynamic model is used to show that the thermodynamic Grüneisen parameter is proportional to the density and independent of temperature. The Mie-Grüneisen equation of state adequately describes the high temperature behavior of molybdenum under the present range of shock loading conditions.

Experimental study on inertial effects in liquid-solid flows

This thesis presents experimental measurements of the rheological behavior of liquid-solid mixtures at moderate Reynolds (defined by the shear rate and particle diameter) and Stokes numbers, ranging from 3 ≤ Re ≤ 1.6 × 10³ and 0.4 ≤ St ≤ 195. The experiments use a specifically designed Couette cylindrical rheometer that allows for probing the transition from transporting a pure liquid to transporting a dense suspension of particles. Measurements of the shear stress are presented for a wide range of particle concentration (10 to 60% in volume) and for particle to fluid density ratio between 1 and 1.05. The effective relative viscosity exhibits a strong dependence on the solid fraction for all density ratios tested. For density ratio of 1 the effective viscosity increases with Stokes number (St) for volume fractions (φ) lower than 40% and becomes constant for higher φ. When the particles are denser than the liquid, the effective viscosity shows a stronger dependance on St. An analysis of the particle resuspension for the case with a density ratio of 1.05 is presented and used to predict the local volume fraction where the shear stress measurements take place. When the local volume fraction is considered, the effective viscosity for settling and no settling particles is consistent, indicating that the effective viscosity is independent of differences in density between the solid and liquid phase. Shear stress measurements of pure fluids (no particles) were performed using the same rheometer, and a deviation from laminar behavior is observed for gap Reynolds numbers above 4× 10³, indicating the presence of hydrodynamic instabilities associated with the rotation of the outer cylinder. The increase on the effective viscosity with Stokes numbers observed for mixtures with φ ≤ 30% appears to be affected by such hydrodynamic instabilities. The effective viscosity for the current experiments is considerably higher than the one reported in non-inertial suspensions.

I. The effect of droplet solidification upon two-phase flow in a rocket nozzle. II. A kinetic theory investigation of some condensation-evaporation phenomena by a moment method

Part I:

The perturbation technique developed by Rannie and Marble is used to study the effect of droplet solidification upon two-phase flow in a rocket nozzle. It is shown that under certain conditions an equilibrium flow exists, where the gas and particle phases have the same velocity and temperature at each section of the nozzle. The flow is divided into three regions: the first region, where the particles are all in the form of liquid droplets; a second region, over which the droplets solidify at constant freezing temperature; and a third region, where the particles are all solid. By a perturbation about the equilibrium flow, a solution is obtained for small particle slip velocities using the Stokes drag law and the corresponding approximation for heat transfer between the particle and gas phases. Singular perturbation procedure is required to handle the problem at points where solidification first starts and where it is complete. The effects of solidification are noticeable.

Part II:

When a liquid surface, in contact with only its pure vapor, is not in the thermodynamic equilibrium with it, a net condensation or evaporation of fluid occurs. This phenomenon is studied from a kinetic theory viewpoint by means of moment method developed by Lees. The evaporation-condensation rate is calculated for a spherical droplet and for a liquid sheet, when the temperatures and pressures are not too far removed from their equilibrium values. The solutions are valid for the whole range of Knudsen numbers from the free molecule to the continuum limit. In the continuum limit, the mass flux rate is proportional to the pressure difference alone.

Advancements in jet turbulence and noise modeling: accurate one-way solutions and empirical evaluation of the nonlinear forcing of wavepackets

Jet noise reduction is an important goal within both commercial and military aviation. Although large-scale numerical simulations are now able to simultaneously compute turbulent jets and their radiated sound, lost-cost, physically-motivated models are needed to guide noise-reduction efforts. A particularly promising modeling approach centers around certain large-scale coherent structures, called wavepackets, that are observed in jets and their radiated sound. The typical approach to modeling wavepackets is to approximate them as linear modal solutions of the Euler or Navier-Stokes equations linearized about the long-time mean of the turbulent flow field. The near-field wavepackets obtained from these models show compelling agreement with those educed from experimental and simulation data for both subsonic and supersonic jets, but the acoustic radiation is severely under-predicted in the subsonic case. This thesis contributes to two aspects of these models. First, two new solution methods are developed that can be used to efficiently compute wavepackets and their acoustic radiation, reducing the computational cost of the model by more than an order of magnitude. The new techniques are spatial integration methods and constitute a well-posed, convergent alternative to the frequently used parabolized stability equations. Using concepts related to well-posed boundary conditions, the methods are formulated for general hyperbolic equations and thus have potential applications in many fields of physics and engineering. Second, the nonlinear and stochastic forcing of wavepackets is investigated with the goal of identifying and characterizing the missing dynamics responsible for the under-prediction of acoustic radiation by linear wavepacket models for subsonic jets. Specifically, we use ensembles of large-eddy-simulation flow and force data along with two data decomposition techniques to educe the actual nonlinear forcing experienced by wavepackets in a Mach 0.9 turbulent jet. Modes with high energy are extracted using proper orthogonal decomposition, while high gain modes are identified using a novel technique called empirical resolvent-mode decomposition. In contrast to the flow and acoustic fields, the forcing field is characterized by a lack of energetic coherent structures. Furthermore, the structures that do exist are largely uncorrelated with the acoustic field. Instead, the forces that most efficiently excite an acoustic response appear to take the form of random turbulent fluctuations, implying that direct feedback from nonlinear interactions amongst wavepackets is not an essential noise source mechanism. This suggests that the essential ingredients of sound generation in high Reynolds number jets are contained within the linearized Navier-Stokes operator rather than in the nonlinear forcing terms, a conclusion that has important implications for jet noise modeling.

A multiple scattering problem

The present work deals with the problem of the interaction of the electromagnetic radiation with a statistical distribution of nonmagnetic dielectric particles immersed in an infinite homogeneous isotropic, non-magnetic medium. The wavelength of the incident radiation can be less, equal or greater than the linear dimension of a particle. The distance between any two particles is several wavelengths. A single particle in the absence of the others is assumed to scatter like a Rayleigh-Gans particle, i.e. interaction between the volume elements (self-interaction) is neglected. The interaction of the particles is taken into account (multiple scattering) and conditions are set up for the case of a lossless medium which guarantee that the multiple scattering contribution is more important than the self-interaction one. These conditions relate the wavelength λ and the linear dimensions of a particle a and of the region occupied by the particles D. It is found that for constant λ/a, D is proportional to λ and that |Δχ|, where Δχ is the difference in the dielectric susceptibilities between particle and medium, has to lie within a certain range.

The total scattering field is obtained as a series the several terms of which represent the corresponding multiple scattering orders. The first term is a single scattering term. The ensemble average of the total scattering intensity is then obtained as a series which does not involve terms due to products between terms of different orders. Thus the waves corresponding to different orders are independent and their Stokes parameters add.

The second and third order intensity terms are explicitly computed. The method used suggests a general approach for computing any order. It is found that in general the first order scattering intensity pattern (or phase function) peaks in the forward direction Θ = 0. The second order tends to smooth out the pattern giving a maximum in the Θ = π/2 direction and minima in the Θ = 0 , Θ = π directions. This ceases to be true if ka (where k = 2π/λ) becomes large (> 20). For large ka the forward direction is further enhanced. Similar features are expected from the higher orders even though the critical value of ka may increase with the order.

The first order polarization of the scattered wave is determined. The ensemble average of the Stokes parameters of the scattered wave is explicitly computed for the second order. A similar method can be applied for any order. It is found that the polarization of the scattered wave depends on the polarization of the incident wave. If the latter is elliptically polarized then the first order scattered wave is elliptically polarized, but in the Θ = π/2 direction is linearly polarized. If the incident wave is circularly polarized the first order scattered wave is elliptically polarized except for the directions Θ = π/2 (linearly polarized) and Θ = 0, π (circularly polarized). The handedness of the Θ = 0 wave is the same as that of the incident whereas the handedness of the Θ = π wave is opposite. If the incident wave is linearly polarized the first order scattered wave is also linearly polarized. The second order makes the total scattered wave to be elliptically polarized for any Θ no matter what the incident wave is. However, the handedness of the total scattered wave is not altered by the second order. Higher orders have similar effects as the second order.

If the medium is lossy the general approach employed for the lossless case is still valid. Only the algebra increases in complexity. It is found that the results of the lossless case are insensitive in the first order of k_imD where k_im = imaginary part of the wave vector k and D a linear characteristic dimension of the region occupied by the particles. Thus moderately extended regions and small losses make (k_imD)² ≪ 1 and the lossy character of the medium does not alter the results of the lossless case. In general the presence of the losses tends to reduce the forward scattering.

Non-Linear Scale Interactions in a Forced Turbulent Boundary Layer

This thesis explores the dynamics of scale interactions in a turbulent boundary layer through a forcing-response type experimental study. An emphasis is placed on the analysis of triadic wavenumber interactions since the governing Navier-Stokes equations for the flow necessitate a direct coupling between triadically consist scales. Two sets of experiments were performed in which deterministic disturbances were introduced into the flow using a spatially-impulsive dynamic wall perturbation. Hotwire anemometry was employed to measure the downstream turbulent velocity and study the flow response to the external forcing. In the first set of experiments, which were based on a recent investigation of dynamic forcing effects in a turbulent boundary layer, a 2D (spanwise constant) spatio-temporal normal mode was excited in the flow; the streamwise length and time scales of the synthetic mode roughly correspond to the very-large-scale-motions (VLSM) found naturally in canonical flows. Correlation studies between the large- and small-scale velocity signals reveal an alteration of the natural phase relations between scales by the synthetic mode. In particular, a strong phase-locking or organizing effect is seen on directly coupled small-scales through triadic interactions. Having characterized the bulk influence of a single energetic mode on the flow dynamics, a second set of experiments aimed at isolating specific triadic interactions was performed. Two distinct 2D large-scale normal modes were excited in the flow, and the response at the corresponding sum and difference wavenumbers was isolated from the turbulent signals. Results from this experiment serve as an unique demonstration of direct non-linear interactions in a fully turbulent wall-bounded flow, and allow for examination of phase relationships involving specific interacting scales. A direct connection is also made to the Navier-Stokes resolvent operator framework developed in recent literature. Results and analysis from the present work offer insights into the dynamical structure of wall turbulence, and have interesting implications for design of practical turbulence manipulation or control strategies.

Fast Lattice Green's Function Methods for Viscous Incompressible Flows on Unbounded Domains

In this thesis, a collection of novel numerical techniques culminating in a fast, parallel method for the direct numerical simulation of incompressible viscous flows around surfaces immersed in unbounded fluid domains is presented. At the core of all these techniques is the use of the fundamental solutions, or lattice Green’s functions, of discrete operators to solve inhomogeneous elliptic difference equations arising in the discretization of the three-dimensional incompressible Navier-Stokes equations on unbounded regular grids. In addition to automatically enforcing the natural free-space boundary conditions, these new lattice Green’s function techniques facilitate the implementation of robust staggered-Cartesian-grid flow solvers with efficient nodal distributions and fast multipole methods. The provable conservation and stability properties of the appropriately combined discretization and solution techniques ensure robust numerical solutions. Numerical experiments on thin vortex rings, low-aspect-ratio flat plates, and spheres are used verify the accuracy, physical fidelity, and computational efficiency of the present formulations.