The growing use of cardiovascular simulations for diagnosis and surgical planning is facilitated by constant advancements in numerical methods and modeling. This dissertation is a collection of three studies, two of which introduced improvements to the current CFD solver that improves the solution accuracy and speed, while the last study uses CFD simulations to perform a parametric optimization study, all in relation to cardiovascular flow. The first study introduced a method to perform cardiovascular CFD simulations in the Fourier domain. The growing use of cardiovascular simulations for diagnosis and surgical planning requires faster computational methods than those currently available. To address this need, we leverage the periodic nature of these flows by discretizing equations in terms of Fourier modes instead of time steps. This approach, known as the harmonic balance method, significantly reduces the size of the discrete problem and hence the simulation cost. In our study, we introduce a harmonic balance finite element solver for simulating physically stable time-periodic flows. The proposed solver is formulated using the baseline Galerkin's method with a least-squares stabilization term. The solver is tested on three physiological cases: a Glenn operation pulmonary arteries flow, a cerebral arteries flow, and a left main coronary arteries flow, against a conventional time-stepping solver. We demonstrate that simulations using the harmonic balance solver all converged around 30 minute, while cases run with the conventional time solver takes more than ten hours, resulting in 10 to 100 times simulation speedup. We compared the solutions between the proposed and conventional solvers and found that the results are very similar, with a margin of around 5%. The second study proposed a new stabilization formulation for the streamline upwind Petrov-Galerkin finite element method. Several finite element methods for simulating incompressible flows rely on the streamline upwind Petrov-Galerkin stabilization (SUPG) term, which is weighted by $\tau_{\mathrm{SUPG}}$. The conventional formulation of $\tau_{\mathrm{SUPG}}$ includes a constant that depends on the time step size, producing an overall method that becomes exceedingly less accurate as the time step size approaches zero. In practice, such method inconsistency introduces significant error in the solution, especially in cardiovascular simulations, where small time step sizes may be required to resolve multiple scales of the blood flow. To overcome this issue, we propose a consistent method that is based on a new definition of $\tau_{\mathrm{SUPG}}$. This method, which can be easily implemented on top of an existing streamline upwind Petrov-Galerkin and pressure stabilizing Petrov-Galerkin method, involves the replacement of the time step size in $\tau_{\mathrm{SUPG}}$ with a physical time scale. This time scale is calculated in a simple operation once every time step for the entire computational domain from the ratio of the L2-norm of the acceleration and the velocity. The proposed method is compared against the conventional method using four cases: a steady pipe flow, a blood flow through vascular anatomy, an external flow over a square obstacle, and a fluid-structure interaction case involving an oscillatory flexible beam. These numerical experiments, which are performed using linear interpolation functions, show that the proposed formulation eliminates the inconsistency issue associated with the conventional formulation in all cases. While the proposed method is slightly more costly than the conventional method, it significantly reduces the error, particularly at small time step sizes. For the pipe flow where an exact solution is available, we show the conventional method can over-predict the pressure drop by a factor of three. This large error is almost completely eliminated by the proposed formulation, dropping to approximately 1% for all time step sizes and Reynolds numbers considered. The third study introduces an algebraic model informed by computational fluid dynamics (CFD) simulations to investigate the performance of the assisted bidirectional Glenn (ABG) operation on a broad range of conditions. The performance of this operation, as measured by the superior vena cava (SVC) pressure, depends on the nozzle area in its ejector pump and the patient's pulmonary vascular resistance (PVR). Using the developed algebraic model to explore this two-dimensional parameter space shows that the ejector pump can create a pressure difference between the pulmonary artery and the SVC as high as 5 mmHg. The lowest SVC pressure is produced at a nozzle area that decreases linearly with the PVR such that, at PVR = 4.2 (Wood units-m2), there is no added benefit in utilizing the ejector pump effect (optimal nozzle area is zero, corresponding to the bidirectional Glenn circulation). At PVR = 2 (Wood units-m2), the SVC pressure can be lowered to less than 4 mmHg by using an optimal nozzle area of approximately 2.5 mm2. Regardless of the PVR, adding a 2 mm2 nozzle to the baseline bidirectional Glenn boosts the oxygen saturation and delivery by at least 15%. The SVC pressure for that 2 mm2 nozzle remains below 14 mmHg for all PVRs less than 7 Wood units-m2. The mechanical efficiency of the optimal designs consistently remains below 30%, indicating the potential for improvement in the future. A good agreement is observed between the algebraic model and high-fidelity CFD simulations.