This paper deals with a multiphysics numerical modelling via finite element method (FEM) of an air-coupled array piezoelectric micromachined ultrasonic transducers (PMUTs). The proposed numerical model is fully 3D with the following features: the presence of the fabrication induced residual stresses, which determine a geometrically non-linear initial deformed configuration of the diaphragms and a remarkable shift of the fundamental frequency; the multiple coupling between different physics, namely electro-mechanical-coupling for the piezo-electric model, acoustic-structure interaction at the acoustic-structure interface and pressure acoustics in the surrounding air. The model takes into account the complete set of PMUTs belonging to the silicon die in a 4 × 4 array configuration and the protective package, as well. The results have been validated by experimental data, in terms of initial static pre-deflected configuration of the diaphragms and frequency response function of the PMUT. The numerical procedure was applied, to analyze different package configurations of the device, to study the influence of the holes on the acoustic transmission in terms of SPL and propagation pattern and consequently extract a set of design guidelines.