The aim of this work has been to implement a range of Large Eddy Simulation (LES) turbulence models into Oasis, a Computational Fluid Dynamic (CFD) solver for incompressible flows based on the Finite Element method (FEM), developed in-house at the University of Oslo. The focus has been on subgrid-scale models that apply the eddy viscosity hypothesis to close the equation set, where both static and dynamic types have been investigated. Verification and assessment of the implementation is performed applying the Method of Manufactured Solutions (MMS), and the classic case of fully developed turbulent channel flow in an x-z periodic channel. The work is further validated through the U.S. Food and Drug Administration (USFDA/FDA)'s computational round robin #1. MMS does in general return good results, where the convergence rate in time is correctly computed to r = 2 for some constructed eddy viscosity expressions, the Smagorinsky model, and the WALE model. As for the channel flow case, the Smagorinsky model returns mean velocity profiles that are closer to those of resolved Direct Numerical Simulations (DNS), compared to what is obtained with under-resolved DNS. The WALE model, the Sigma model, and the Dynamic Smagorinsky model, which all have the correct wall behaviour, return results where little or no improvements are seen compared to the under-resolved profiles. It turns out that a net contribution of eddy viscosity close to the wall is extremely important for this specific case, as this in some sense controls the whole quality of the simulation. The Sigma model and the Dynamic Smagorinsky model do, in despite of small improvements to the mean velocities, return good profiles for some selected Reynolds stresses. For the FDA case good results are obtained for both uniform and non-uniform meshes for all LES models, where the mean velocity profiles in general are between 50 to 80 percent closer to experimental data, compared to what is obtained with under-resolved DNS. The only model that, as expected, erroneously predicts the breaking position of the jet is the Smagorinsky model. The best validation metric is obtained with the WALE model on the non-uniform mesh of three million cells; on the other hand, the best percentage improvement over under-resolved DNS is obtained for the Dynamic Smagorinsky model with Lagrangian averaging on the non-uniform mesh of two million cells.