Studies on the haemodynamics of human circulation are clinically and scientifically important. In order to investigate the effect of deformation and aggregation of red blood cells (RBCs) in blood flow, a computational technique has been developed by coupling the interaction between the fluid and the deformable RBCs. Parallelization was carried out for the coupled code and a high speedup was achieved based on a spatial decomposition. In order to verify the code's capability of simulating RBC deformation and transport, simulations were carried out for a spherical capsule in a microchannel and multiple RBC transport in a Poiseuille flow. RBC transport in a confined tube was also carried out to simulate the peristaltic effects of microvessels. Relatively large-scale simulations were carried out of the motion of 49,512 RBCs in shear flows, which yielded a hematocrit of 45%. The large-scale feature of the simulation has enabled a macroscale verification and investigation of the overall characteristics of RBC aggregations to be carried out. The results are in excellent agreement with experimental studies and, more specifically, both the experimental and simulation results show uniform RBC distributions under high shear rates (60-100/s) whereas large aggregations were observed under a lower shear rate of 10/s.