Immersed boundary (IB) methods are attractive due to their ability to simulate flow over complex geometries on a simple Cartesian mesh. Unlike conformal grid formulation, the mesh does not need to conform to the shape and orientation of the boundary. This eliminates the need for complex mesh and/or re-meshing in simulations with moving/morphing boundaries, which can be cumbersome and computationally expensive. However, the imposition of boundary conditions in IB methods is not straightforward and numerous modifications and refinements have been proposed and a number of variants of this approach now exist. In a nutshell, IB methods in the literature often suffer from numerical oscillations, implementation complexity, time-step restriction, burred interface, and lack of generality. This limits their ability to mimic conformal grid results and enforce Neumann boundary conditions. In addition, there is no generic IB capable of solving flow with multiple potentials, closely/loosely packed structures as well as IBs of infinitesimal thickness. This dissertation describes a novel 2$ ^{\text{nd}} $ order direct forcing immersed boundary method designed for simulation of two- and three-dimensional incompressible flow problems with complex immersed boundaries. In this formulation, each cell cut by the IB is reshaped to conform to the shape of the IB. IBs are modeled as a series of 2D planes in 3D space that connect seamlessly at the edges of the cut cells, in a way that mimics conformal grid. IBs are represented in a continuous and consistent fashion from one cell to another, thus eliminating spatial pressure oscillations originating from inconsistent description of the IB as well as the traditional stair-step problem, leading to a more accurate resolution of the boundary layer. Boundary conditions are enforced at the exact location of the IB devoid of interpolation, which guarantees sound simulations even on grids with high aspect ratio, and enables simulations of flow packed with multiple IBs in close proximity. Boundary conditions for each phase across the IB are enforced independently, yielding a unique capability to solve flows with zero-thickness IBs. Simulations of a large number of 2D and 3D test cases confirm the prowess of the devised immersed boundary method in solving flows over multiple loosely/closely-packed IBs; stationary, moving and highly morphing IBs; as well as IBs with zero-thickness. Extension of the proposed scheme to solve flow with multiple potentials is demonstrated by simulating transfer and transport of a passive scalar from an array of side-by-side and tandem cylinders in cross-flow. Aquatic vegetation represented by a colony of circular cylinders with low to high solid fraction is simulated to showcase the prowess of the current numerical technique in solving flow with closely packed structures. Aquatic vegetation studies are extended to a colony of flat plates with different orientations to show the capability of the developed method in modeling zero-thickness structures.
