Dense suspensions can be found in various industrial and natural processes. A relatively new technique uses the principle of additive manufacturing to produce products from a wide variety of materials by printing with dense suspensions. To reach a high print quality the suspension rheology must be understood very well. Current knowledge about suspension rheology however lacks the capability of predicting the exact behaviour of a predefined suspension, especially for dense suspensions. Direct numerical simulation (DNS) in combination with a second order accurate immersed boundary method (IBM) (Breugem, 2012) can be a useful tool to research suspension rheology. However, currently no numerical results are known for suspensions close to the jamming limit produced with this method. The current work validates the capability of the IBM to simulate dense suspensions by comparing produced results with existing numerical and experimental data. This is done by simulating a plane Couette flow for a range of particle volume fractions 𝜙=0.2−0.6, all of which are simulated with two friction coefficients (𝜇𝑐=0 and 0.39). Furthermore, the focus is on Stokes flow of neutrally buoyant non-colloidal suspensions with monodisperse spherical particles and the channel height was chosen equal to 13.5 particle diameters. DNS has been used to analyse the suspension rheology in terms of mean concentration profiles, velocity profiles, interactions in the microstructure and particle stress profiles. Steady-state concentration profiles of the simulated cases show a particle layering effect close to the confining walls. This layering effect alters the suspension rheology significantly and for that reason the wall regions are analysed separately from the core region. For both regions the microstructure, the relative viscosity and the normal particle stresses are analysed. The results agree well with existing numerical results (Gallier et al., 2016) (Yeo & Maxey, 2010a). Comparison with experimental work (Dbouk et al., 2013) (Zarraga et al., 2000) shows that results for the relative viscosities of the unlayered core regions are lower in general, reasons for this difference can be higher friction factors or the use of non-spherical particles in experiments. However, the same asymptotic trend is observed for increasing 𝜙. The maximum packing fraction 𝜙𝑚 was found by fitting the relative viscosities to the Marron & Pierce equation (Maron & Pierce, 1956), which gave 𝜙𝑚=0.69 for 𝜇𝑐=0 and 𝜙𝑚=0.635 for 𝜇𝑐=0.39. These results for 𝜙𝑚 are in good agreement with results from (Gallier et al., 2014). In general the IBM turns out to be capable of reproducing existing numerical data. Furthermore, results have been obtained for suspensions closer to the jamming limit than known so far in numerical work on this particular flow regime. Besides that, the present results are obtained with an advanced soft-sphere collision model, including lubrication corrections for close approach of particles, which has been extensively validated with collision experiments in a previous study. Comparison of the results from this work with experimental data shows larger differences. The reason for this can be that the suspensions in experimental set-ups deviate from the idealized suspension in the simulations. Besides that experimental results differ significantly from each other, indicating that differences also exist between the experimental suspensions. These differences have to be clearly defined in order to make a valuable comparison. Therefore it is currently difficult to determine how accurately the IBM simulates suspension rheology.