Abstract: In this paper, we propose the modeling and simulation of a flapping bird based on the immersion boundary method. Firstly, the accuracy of the numerical method is tested. And then the heave motion and pitch motion of the NACA 0012 airfoil are studied. Finally, the effect of asymmetry flapping time of a three-dimensional flapping wing is studied. The results show that the simulation of flapping airfoil by immersion boundary method is in good agreement with available results in the literature. The propulsive capacity of heave motion is determined by the size and position of the leading edge vorticity of the airfoil. The peak value of the propulsive efficiency of heave motion is mainly concentrated when 0.3≤St≤0.4. When the heave motion is coupled with the pitch motion, the propulsive efficiency reaches its peak when the pitch angle is 25° and the phase difference is 85°. The results from 3D simulations reveal that appropriately increasing the downbeat speed of the wings provides more lift and consumes more energy. The research results can provide reference for flutter parameter setting of micro flapping wing aircraft.