We present a wavelet-based adaptive method for computing 3D multiscale flows in complex, time-dependent geometries, implemented on massively parallel computers. While our focus is on simulations of flapping insects, it can be used for other flow problems. We model the incompressible fluid with an artificial compressibility approach in order to avoid solving elliptical problems. No-slip and in/outflow boundary conditions are imposed using volume penalization. The governing equations are discretized on a locally uniform Cartesian grid with centered finite differences, and integrated in time with a Runge–Kutta scheme, both of 4th order. The domain is partitioned into cubic blocks with different resolution and, for each block, biorthogonal interpolating wavelets are used as refinement indicators and prediction operators. Thresholding the wavelet coefficients allows to generate dynamically evolving grids, and an adaption strategy tracks the solution in both space and scale. Blocks are distributed among MPI processes and the grid topology is encoded using a tree-like data structure. Analyzing the different physical and numerical parameters allows us to balance their errors and thus ensures optimal convergence while minimizing computational effort. Different validation tests score accuracy and performance of our new open source code, WABBIT. Flow simulations of flapping insects demonstrate its applicability to complex, bio-inspired problems.