This paper researched into the harmonic and anharmonic underwater flapping foil propulsion systems to improve the efficiency of these bioinspired propulsors. The angle of attack, the pitching angle, the heaving amplitude, and the phase difference are parametrically investigated in this paper. A rigid two-dimensional NACA (National Advisory Committee for Aeronautics) 0012 airfoil is modeled with the aid of a commercial computational fluid dynamics software, FINE™/Marine. Unsteady Reynolds Average Navier-Stokes (URANS) equation is solved together with dynamic mesh to simulate the foil motion. The investigation first verifies the reliability of the developed modeling method against the benchmark data. Then, the systematic investigation is conducted and identifies that the heaving amplitude is most influential factor for the propulsion efficiency. Secondly, phase difference also has a significant influence on efficiency, but this effect is related to the reference working condition, which needs further study. Then, the pitching amplitude has little effect on the maximum efficiency value of flapping foil, while it will affect its optimal speed range. When the heaving amplitude ratio reaches 3 and the corresponding maximum angle of attack is about 9°, the maximum efficiency can reach 87%. The effect of anharmonic motion on the efficiency is very small and varies with the St number, but in summary, it can maintain the peak efficiency over a wider range of operations. In addition, the force and flow field characteristics of different efficiency points are compared and analyzed to distinguish their corresponding relationship with the propulsion efficiency.