In multiphase flows, accurately modeling the interaction between the liquid phase of complex fluids and a porous medium of solid spheres poses a fundamental challenge. The dynamics of moderately dense non-colloidal suspensions constituted by static random arrays of mono-disperse spherical particles in non-linear viscoelastic fluids is studied numerically. This numerical study consists of about 9000 different systems, in which the volume fraction ϕ (0.04≤ϕ≤0.2) of the dispersed solid phase, the Reynolds number Re(5≤Re≤50), the solvent viscosity ratio β(0.05≤β≤0.9), the Weissenberg number Wi(0.5≤Wi≤4), and the mobility parameter of the Giesekus model α (0.1≤α≤0.5) were varied to understand the particle's interactions with the viscoelastic suspending fluid. We aim to investigate the relationship between the volume fraction of the dispersed solid phase and the non-linear rheology of shear-thinning viscoelastic fluids with the normalized average drag force ⟨F⟩. In addition, by assessing the flow patterns predicted numerically, we were able to provide a characterization of the velocity and stress fields as a function of the simulation parameters.