Blood-based microfluidic devices typically have characteristic dimensions ranging from tens of micrometers to several hundreds. At these scales, blood exhibits strong non-Newtonian behaviour, such as shear-thinning, as well as multiphase dynamics, such as the migration of Red Blood Cells (RBCs) toward the centreline of microchannels and the formation of a cell-free layer. The development process of such microfluidic devices is mainly driven by performance metrics such as separation efficiency and hematocrit distribution. It is thus crucial that numerical methods accurately capture these flow behaviours to predict device performance and inform the design process. This work presents the numerical investigation of blood flow in microfluidic devices using a modified Eulerian-Eulerian two-phase model within the OpenFOAM framework (twoPhaseEulerFoam). Blood is assumed to be a suspension of RBCs within plasma, where both phases are considered as interpenetrating continua with their own conservation equations. Specifically, blood plasma is treated as a Newtonian fluid, while the RBCs are modelled as a shear-thinning fluid with a hematocrit-based viscosity model. The model is applied to benchmark microfluidic junction geometries (T- and Y-) typically existing in plasma separation devices. Its accuracy is then evaluated by comparing the simulated results against numerical and experimental data available in the literature.