Many scientific and technical problems which concern the dynamics of complex fluids such as multi-phase-flow and realistic flow in porous and granular media deal with the interaction between fluids and particles, rather than with the dynamics of the fluid alone. The research of how the surrounding fluid affects the dynamics of particles, or how to deal with the problem computationally for the microscopic level is still at the beginning. The aim of this study is to develop a microscopic simulation method (fluid goes around the particles) where granular particles can be simulated inside fluids to study those problems. This is done by combining the simulation method for granular particles with the simulation method for the incompressible Newtonian fluid. The granular particles are implemented via the discrete element method (DEM) where the elastic contact force between two undeformed contacting polygonal particles is proportional to the overlap area ("hard particle, soft contact"). The Gear Predictor-Corrector of 2nd-order (BDF2) is used as the time integrator to solve the equations of motion of the particles. For the fluid phase, the implementation of the incompressible Navier-Stokes equations via the Galerkin finite element method (FEM) is formulated as differential algebraic equations (DAE) with the pressures as the Lagrange parameters. The time integration is again via the BDF2 while the resulting non-linear equations are solved via the Newton-Raphson methods. The spatial discretization is via the Taylor-Hood elements from Delaunay triangulations with additional post-processing with the relaxation algorithm. The coupling of the DEM for the granular particles and the FEM for the fluid is via appropriate boundary conditions and the drag force (computed by the integration of the fluid stress tensor over the particle's surface). This is being verified via the computation of wall correction factors of a sinking particle. The fluid simulation is extended to a simulation of free surfaces where the motion of the surface is integrated out according to the velocity on the surface which is obtained from the FEM-scheme. The second-order Adams-Bashforth method turns out to be the most suitable integrator for the surface motion. Compared to conventional efforts, which try to solve partial differential equations for the motion of the surface, the additional effort in our method with respect to new data structures etc. is minimal. The free surfaces code is verified by simulating the collapse of a water column. For the speed of the wavefronts, excellent agreement is obtained for large viscosity with the lubrication approximation. The agreement of the results with the experimental data for water is a further gratifying result. Two numerical experiments are conducted using the DEM-FEM code: one with a rather slow dynamics, another one relatively more "violent". The compaction simulation has shown that the addition of fluid to a granular assembly can increase the sound velocity in the system, compared to the dry case. The high viscosity slowed down the compaction, irrespective whether the system was tapped only on the ground or on the whole boundary. The granular column simulations show that for systems immersed under fluids, rolling of particles becomes less important than for the corresponding dry systems.