In this study, a three-dimensional artificial compressibility solver based on the average-state Harten-Lax-van Leer-Contact (HLLC) [13] type Riemann solution is first proposed and developed to solve the time-dependent incompressible flow equations. To implement unsteady flow calculations, a dual time stepping strategy including the LU decomposition method is used in the pseudo-time iteration and the second-order accurate backward difference is adopted to discretize the unsteady flow term. Also a third-order accurate HLLC numerical flux is derived for approximating the inviscid terms. To verify numerical accuracy, flows over a lid-driven cavity and an oscillating flat plate are chosen as the benchmark tests. In addition, the current solver is extended to solve blood flows in a realistic human aorta measured from MRI (Magnetic Resonance Imaging). The simulation geometry was derived from a three-dimensional reconstruction of a series of two-dimensional slices obtained in vivo. Numerical results demonstrate wall stresses were highly dynamic, but were generally high along the outer wall in the vicinity of the branches and low along the inner wall, particularly in the descending aorta. The maximum wall stress distribution is presented on the aortic arch in the systole. In addition, extensive counter-clockwise secondary flows and three-dimensional helical vortex influenced considerably by the presence of vessel contraction, torsion and the branches were shown in the descending aorta in the late systole and early diastolic cycles.