A Navier-Stokes code based on an unstructured finite volume method is used to simulate the MIT flapping foil experiment. A low Reynolds number ${kappa}-{varepsilon}$ turbulence model is used to close the Reynolds averaged Navier-Stokes equations. Computations are carried out for the whole experimental domain involving two flapping foils and a downstream hydrofoil. The computational domain is meshed with unstructured quadrilateral elements, partly structured. Numerical solutions show good agreement with experiment. The first harmonics of the velocity in the boundary layer shows local peak value inside the boundary layer and also local minimum near the edge of boundary layer. It is intensified as it develops along the blade surface. This is shown to be caused as the unsteadiness inside the boundary layer is being convected at a speed less than the free stream value. It is also shown that there is negligible mixing of the unsteadiness between the boundary layer and the free stream