Determination of the relation between the bulk or rheological properties of a particle suspension and its microscopic structure is an old and important problem in physical science. In general, the rheology of particle suspension is quite complex, and the problem becomes even more complicated if the suspending particle is deformable. Despite these difficulties, a large number of theoretical and experimental investigations have been devoted to the analysis and prediction of the rheological behavior of particle suspensions. However, among these studies there are very few investigations that focus on the role of particle deformability. A novel method for full coupling of the fluid-solid phases with sub-grid accuracy for the solid phase is developed. In this method, the flow is computed on a fixed regular 'lattice' using the lattice Boltzmann method (LBM), where each solid particle, or fiber, is mapped onto a Lagrangian frame moving continuously through the domain. The motion and orientation of the particle are obtained from Newtonian dynamics equations. The deformable particle is modeled by the lattice-spring model (LSM).The fiber deformation is calculated by an efficient flexible fiber model. The no-slip boundary condition at the fluid-solid interface is based on the external boundary force (EBF) method. This method is validated by comparing with known experimental and theoretical results. The fiber simulation results show that the rheological properties of flexible fiber suspension are highly dependent on the microstructural characteristics of the suspension. It is shown that fiber stiffness (bending ratio BR) has strong impact on the suspension rheology in the range BR < 3. The relative viscosity of the fiber suspension under shear increases significantly as BR decreases. Direct numerical simulation of flexible fiber suspension allows computation of the primary normal stress difference as a function of BR. These results show that the primary normal stress difference has a minimum value at BR ∼ 1. The primary normal stress differences for slightly deformable fibers reaches a minimum and increases significantly as BR decreases below 1. The results are explained based on the Batchelor's relation for non-Brownian suspensions. The influence of fiber stiffness on the fiber orientation distribution and orbit constant is the major contributor to the variation in rheological properties. A least-squares curve-fitting relation for the relative viscosity is obtained for flexible fiber suspension. This relation can be used to predict the relative viscosity of flexible fiber suspension based on the result of rigid fiber suspension. The unique capability of the LBM-EBF method for sub-grid resolution and multiscale analysis of particle suspension is applied to the challenging problem of platelet motion in blood flow. By computing the stress distribution over the platelet, the "blood damage index" is computed and compared with experiments in channels with various geometries [43]. In platelet simulation, the effect of 3D channel geometry on the platelet activation and aggregation is modeled by using LBM-EBF method. Comparison of our simulations with Fallon's experiments [43] shows a similar pattern, and shows that Dumont's BDI model [40] is more appropriate for blood damage investigation. It has been shown that channels with sharp transition geometry will have larger recirculation areas with high BDI values. By investigating the effect of hinge area geometry on BDI value, we intend to use this multiscale computational method to optimize the design of Bileaflet mechanical heart valves. Both fiber simulations and platelet simulations have shown that the novel LBM-EBF method is more efficient and stable compare to the conventional numerical methods. The new EBF method is a two-Cway coupling method with sub-grid accuracy which makes the platelet simulations possible. The LBM-EBF is the only method to date, to the best of author's knowledge, that can simulate suspensions with large number of deformable particles under complex flow conditions. It is hoped that future researchers may benefit from this new method and the algorithms developed here.