This thesis considers three aspects of the numerical simulations, which are

coupling, conservation, and performance. We conduct a project and address

one challenge from each of these aspects.

We propose a novel penalty force to enforce contacts with accurate Coulomb

friction. The force is compatible with fully-implicit time integration and the

use of optimization-based integration. In addition to processing collisions

between deformable objects, the force can be used to couple rigid bodies to

deformable objects or the material point method. The force naturally leads to

stable stacking without drift over time, even when solvers are not run to

convergence. The force leads to an asymmetrical system, and we provide a

practical solution for handling these.

Next we present a new technique for transferring momentum and velocity between

particles and MAC grids based on the Affine-Particle-In-Cell (APIC) framework

previously developed for co-located

grids. We extend the original APIC paper and show that

the proposed transfers preserve linear and angular momentum and also satisfy

all of the original APIC properties.

Early indications in the original APIC paper suggested that APIC might be

suitable for simulating high Reynolds fluids due to favorable retention of

vortices, but these properties were not studied further. We use two

dimensional Fourier analysis to investigate dissipation in the limit $\dt=0$.

We investigate dissipation and vortex retention numerically to quantify the

effectiveness of APIC compared with other transfer algorithms.

Finally we present an efficient solver for problems typically seen in

microfluidic applications.

Microfluidic ``lab on a chip'' devices are small devices that operate on small

length scales on small volumes of fluid. Designs for microfluidic chips are

generally composed of standardized and often repeated components connected by

long, thin, straight fluid channels. We propose a novel discretization

algorithm for simulating the Stokes equations on geometry with these features,

which produces sparse linear systems with many repeated matrix blocks. The

discretization is formally third order accurate for velocity and second order

accurate for pressure in the $L^\infty$ norm. We also propose a novel linear

system solver based on cyclic reduction, reordered sparse Gaussian elimination,

and operation caching that is designed to efficiently solve systems with

repeated matrix blocks.

When multiple fluids in a single system interact, they can mix and separate according to their physical properties. We can describe a two-fluid system of this kind using the phase-field model. This model uses a scalar field to represent the volume concentration or concentration difference of one of the fluids. By varying this field we can describe phase mixtures, and the interface between phases is represented by the field smoothly transitioning between extreme values. The Cahn-Hilliard equation is a fourth-order partial differential equation that can be used to describe the change in this phase-field through separation. This method acts to minimize an energy potential function with minima where the mixture is fully separated. We couple the Cahn-Hilliard equation to the classic Navier-Stokes equations to include fluid flow. We develop two novel discretization methods for this system. The first method is a simplified time integration scheme based on splitting the fourth-order equation into a pair of second-order Helmholtz equations that can be solved as a pair of linear systems. We propose a method of coupling this discretization to an incompressible Navier-Stokes equation that maintains consistency between mass and momentum. We also propose a novel surface tension force discretization that limits surface forces to the interface region. This surface tension can be controlled with a parameter that allows for transitioning between a sharp interface formulation and a smooth continuum surface force. A drawback of this scheme is that it relies on a mass redistribution post-processing step to maintain the phase variable in the physical bounds. The second method is developed to intrinsically maintain the phase variable inside a physically consistent range instead of requiring a post-processing step. We formulate a piecewise energy potential function with energy barriers at both extremes. We combine this with a nonlinear second-order solve that is implicit in the energy potential. We show how to formulate this nonlinear equation as an optimization problem to be solved with a damped Newton's method. We compare both methods to previous work using a variety of experiments, and demonstrate that they are second-order accurate.