Gravitational waves are one of the most exciting predictions of General Relativity. Due to their compactness, binary black hole (BBH) systems are likely sources for gravitational waves detectable by ground-based interferometric detectors such as Advanced LIGO. Distinguishing a true signal from noise requires accurate models of gravitational waves from BBH systems. In the last decade, numerical relativity has been instrumental in generating these models. The Spectral Einstein Code (SpEC) developed by the SXS Collaboration has been used to investigate many aspects of BBH systems. Generically, the spins of the black holes in the binary are misaligned with the orbital angular momentum, causing the orbital plane and the spins to precess. Such systems are of particular importance as they exhibit interesting dynamics and large modulations of the gravitational waveform.

In this thesis, we explore the problem of modelling precessing BBH by extending the capabilities of SpEC both in the construction of initial data and in dynamical evolution. We thereafter compare numerical relativity results to Post-Newtonian theory to assess the accuracy of precessing dynamics in Post-Newtonian theory.

We begin this dissertation by examining the problem of robustly constructing initial data for high mass ratio, high spin precessing BBH. We discuss many technical improvements that now enable the construction of constraint-satisfying initial data for a much larger region of the parameter space.

Next, we discuss the implementation of an important technical improvement that permits the evolution of arbitrarily precessing BBH systems; in particular, those where the orientation of the orbital plane changes by more than 90 degrees.

We then compare the results of precessing BBH simulations done with SpEC to PN theory. We find generally good agreement between PN and NR precession dynamics, supporting the creation of phenomenological waveforms constructed using rotated nonprecessing waveforms.

Finally, we examine a case of transitional precession in NR, where the total angular momentum changes drastically during the inspiral, and once more find good agreement with PN predictions.