Flapping foils located beneath or to the side of the hull of the ship can be used as unsteady thrusters, augmenting ship propulsion in waves. The basic setup is composed of a horizontal wing, which undergoes an induced vertical motion due to the ship’s responses in waves, while the self-pitching motion of the wing is controlled. Flapping foil thrusters can achieve high level of thrust as indicated by measurements and numerical simulations. Due to the relatively small submergence of the above biomimetic ship thrusters, the free-surface effects become significant. In the present work, the effect of the free surface on the performance of flapping foil thruster is assessed by means of two in-house developed computational models. On one hand, a cost-effective time-domain boundary element method (BEM) solver exploiting parallel programming techniques and general purpose programming on graphics processing units (GPGPU) is employed, while on the other hand a higher fidelity RANSE finite volume solver implemented for high performance computing (HPC) is used, and comparative results are presented. BEM and RANSE calculations present quite similar trends with respect to mean submergence depth, presenting 12%, 28%, and 18% of differences concerning the mean values of lift, thrust, and moment coefficients, respectively. The latter differences become very small after enhancement of the BEM model to include viscous corrections. Useful information and data are derived supporting the design of the considered biomimetic thrusters, for moderate submergence depths and conditions characterized by minor flow separation effects.