After billions of years of evolution, many creatures employing flapping wing in nature tend to have excellent flight capabilities. To understand the bionic wings flow mechanism will be helpful to design high performance underwater vehicles and new conception aircrafts. The geometric parameters, kinematic parameters and flow parameters have great effects on the bionic wings thrust performance. Facing the diverse parameters, it's very difficult to explore the three-dimensional (3D) bionic flapping wing flow mechanism with traditional numerical simulation method. In this paper, a general large-scale parallel solver using Immersed Boundary-Lattice Boltzmann Method (IB-LBM) was developed. The evolution procedures of the 3D flapping wing leading edge vortex and wake flow vortex structures were analyzed in detail. Our study explained the 3D flapping wing thrust performance variation with different wing shapes, aspect ratios and pitch-bias angles of attack. Using Chinese supercomputer TianHe-II presents a wide range of possibilities for the further development of parallel IB-LBM, employing tens of millions grids will help us to obtain more complete and accurate 3D flapping wing flow field information. It's indicated that the obtuse wing has the best thrust performance compared with other sharp wing shapes. With the increase of the aspect ratio, the thrust coefficient of flapping wing increases firstly and then decreases, and with pitch-bias angles of attack increases, the thrust coefficient decreases quickly or even shown resistance phenomenon at the large pitch-bias angel of attack. The discussion of these parameters will provide a theoretical basis for improving flapping-like vehicles propulsive performance.