This study is focused on performing simulations on a NACA 0015 flapping airfoil undergoing heaving and pitching motion used for energy harvesting. The conditions studied are reduced frequency of k = fc/U_∞ = 0.2, pitching amplitude of θ₀ = 70° and heaving amplitude of h₀/c = 1, where f is an oscillating frequency, c is the chord length, U_∞ is the freestream velocity. Results are compared using different turbulence models, the Spalart-Allmaras (SA), k − ω SST, Reynolds Stress, and Scale Adaptive Simulation (SAS) models. Key differences in performance that are evaluated include lift, drag, vortex formation, and computational time among the models for Reynolds' numbers ranging from 7,000 to 350,000. Simulations indicate that near the maxima heaving positions when the pitching motion begins to change direction, there is a peak lift coefficient which indicates large temporal oscillations. These oscillations are linked to vortices shed from the leading and trailing edges of the airfoil. However, differences in the details of these oscillations occur for each model for all Reynolds numbers. Higher peak values of lift force coefficient are obtained for the SAS model, which are up to 7% higher relative to the SA model. These differences lead to different predictions of the time-averaged power coefficient. With identical mesh and time step conditions, the computational time for SA and k − ω SST is almost identical, while the computational time for SAS model was 10x longer.