Complex fluids comprise two phases of materials: a solvent phase and a non-continuum phase such as microscopic particles or polymers. Suspensions of colloids -- microscopic particles small enough to undergo Brownian motion -- are an important example and model system for understanding general complex fluids. The particles are distributed in the suspending solvent, forming a "microstructure". The formulation of detailed models for the dynamics of condensed soft matter including colloidal suspensions and other complex fluids requires accurate description of the physical forces between microstructural constituents. In dilute suspensions, pair-level interactions are sufficient to capture hydrodynamic, interparticle, and thermodynamic forces. In dense suspensions, many-body interactions must be considered. Prior analytical approaches to capturing such interactions such as mean-field approaches replace detailed interactions with averaged approximations. However, long-range coupling and effects of concentration on local structure, which may play an important role in e.g. phase transitions, are smeared out in such approaches. An alternative to such approximations is the detailed modeling of hydrodynamic interactions utilizing precise couplings between moments of the hydrodynamic traction on a suspended particle and the motion of that or other suspended particles. For two isolated spheres, a set of these functions was calculated by Jeffrey and Onishi, Kim and Karrila, and Jeffrey. Along with pioneering work by Batchelor and Green, these are the touchstone for low-Reynolds-number hydrodynamic interactions and have been applied directly in the solution of many important problems related to the dynamics of dilute colloidal dispersions. The Reynolds number, Re, is the dimensionless strength of flow inertia relative to viscous stress, $Re = \rho U L / \eta$, where $\rho$ is the fluid density, $\eta$ is the fluid viscosity, $U$ is the characteristic speed of the flow, and $L$ is the characteristic length scale such as channel width, particle size, or a surface dimension. The overall goal of this dissertation is development of new theoretical and computational models for many-body interactions in concentrated colloidal suspensions. In the first part of the dissertation, we extend the theoretical framework of the pair-level mobility functions to concentrated systems, utilizing a new stochastic sampling technique we developed in order to rapidly calculate an analogous set of mobility functions describing the hydrodynamic interactions between two hard spheres immersed in a suspension of arbitrary concentration. This was carried out utilizing Accelerated Stokesian Dynamics simulations. These mobility functions provide precise, radially dependent couplings of hydrodynamic traction moments to particle velocity derivatives, for arbitrary colloid volume fraction $\phi=4\pi\eta a^3/n$, where $\eta$ is the suspending solvent viscosity, $a$ is the hydrodynamic particle radius, and $n=N/V$ is the number density of colloids in the volume $V$ of solvent. The most significant outcome of this work was to show that hydrodynamic entrainment of one particle by the disturbance flow created by a forced particle decays algebraically just as slowly as a dilute suspension: it decays linearly in the inverse separation distance, ~1/r, between the forced particle and a test particle a distance r away. In contrast with prior assertions in the literature that crowding screens hydrodynamic interactions, this finding holds for volume fractions $0.05\leq\phi\leq0.5$. At these higher concentrations, the coefficients also reveal liquid-like structural effects on pair mobility at close separations. These results confirm that long-range many-body hydrodynamic interactions are an essential part of the dynamics of concentrated systems and that care must be taken when applying renormalization schemes. In Chapter 3, this stochastic technique developed in Chapter 2 is utilized to compute pair mobility couplings between stresslet and straining motion, higher-order traction moments required for many-body couplings away from equililbrium. Thus, the couplings presented in these two chapters constitute a set of orthogonal coupling functions that are utilized to compute equilibrium properties in suspensions at arbitrary concentration and to solve many-body hydrodynamic interactions analytically. We utilize these concentrated mobility functions to extend recently developed dilute theory for the stress in concentrated colloidal suspensions, and compare those results to direct measurement in dynamic simulation. Particle-phase stress in flowing colloidal suspensions represents the extent to which energy is stored entropically or enthalpically by microstructural distortion, minus the energy dissipated by viscous drag. It was recently shown that pair-level hydrodynamic interactions suppress energy storage in dilute dispersions of repulsive hard spheres, with corresponding changes in normal stresses, normal stress differences, and osmotic pressure. In the pair limit, particle roughness or Brownian motion leads to non-Newtonian rheology, whereas pair-level hydrodynamic interactions preserve the structural fore-aft symmetry and contributes to Newtonian rheology. However, in concentrated suspensions where three-body hydrodynamic interactions matter, how three-body interactions mechanistically change rheology is not fully understood. In Chapter 4, we investigate the dependence of non-Newtonian rheology of colloidal dispersions on particle concentration, with a focus on the role of played by pair-level and three-body particle interactions. To do so, we utilize Accelerated Stokesian Dynamics to simulate the evolution of particle-phase stress in concentrated colloidal dispersions undergoing microrheological flow, obtaining detailed measurements of particle structure and dynamics as they evolve with particle concentration, and connect these to changes in non-Newtonian rheology. We find that, in contrast to dilute suspensions where pair hydrodynamic interactions suppress normal stresses and osmotic pressure, suspension stress is enhanced by three-body hydrodynamic interactions, where their role in promoting structural asymmetry produces a concentration-dependent non-Newtonian rheology. The loss of fore-aft symmetry of a pair trajectory in the presence of a third particle inspires us to seek a concentration-dependent hydrodynamic coupling that predicts non-Newtonian behavior under the influence of many-body interactions. To do so, we utilize the set of concentrated pair hydrodynamic functions developed in Chpater 2 to extend dilute theory to concentrated suspension via scaling arguments. Scaling theory and simulations show excellent agreement for the normal stress difference and osmotic pressure, providing support for the idea that, in the presence of a third particle, the transverse displacement of a pair encounter is responsible for the non-Newtonian rheology. We also find that microviscosity is enhanced as concentration increases, and we developed a scaling theory which collapses data of different volume fractions onto the dilute theory. Many-body hydrodynamic interactions play a crucial role in colloidal suspensions , but are notoriously difficult and expensive to model computationally. The Stokesian Dynamics algorithm is one approach to involve many-body hydrodynamic interactions that couple fluid and particle motion, with the primary advantage that detailed fluid motion is not explicitly computed, saving considerable computational expense. However, since hydrodynamic disturbance flows propagate between the particles in an infinite hierarchy of reflections, modeling these so-called many-body hydrodynamic interactions requires multiple matrix operations that drive the primary computational expense of the algorithm. Separation of interactions into near-field and far-field calculations permits analytical treatment of the former and faster computation overall. But in its most optimized form -- Accelerated Stokesian Dynamics (ASD), the algorithm is serial and can handle at most a few thousand particles. However, ever-increasing interest in understanding properties of large-scale systems in colloidal suspensions, such as colloidal gels, has demanded development of techniques to simulate $O(10^5)$ particles. Currently techniques to handle such large-scale systems, such as LAMMPS, neglect hydrodynamic interactions. Other techniques which consider hydrodynamic interactions between colloids utilizing GPU architecture are limited to pair-level hydrodynamic interactions. However, many open questions in the complex fluids literature -- such as the mechanism of the colloidal glass transition, the role of hydrodynamics in colloidal gel collapse, and more -- require accurate modeling of concentrated, hydrodynamically interacting colloids. These systems require $O(10^5)$ to $O(10^6)$ particles to evolve a statistically representative set of network structures. In chapter 5, we present a scalable and parallel algorithms of Accelerated Stokesian Dynamics in a distributed memory architecture to fulfill the requirement of simulating large scale particle systems suspended in Stokes flow. We solve the ASD equations with Krylov subspace methods for the sparse near-field two-body interactions and for the full many-body far-field interactions. The far-field action is matrix-free and is based on Fast Fourier Transforms. We treat in details the cases of shear flow, brownian motion, and the introduction of attractive inter-particle force. The parallelization of the different phases of the algorithm are presented and analyzed and we show scaling up to 8192 processors for 819,200 particles.