Fractured-vuggy carbonate reservoirs are usually complex in structure, containing not only matrix and fractures but also irregularly shaped vugs ranging in diameter from millimeters to meters. These characteristics lead to complex fluid flow in such media, with significant deformation of fractures and vugs. The conventional multiphase flow theory based on Darcy's law can't accurately simulate the multiphase flow in such multiscale media. In this paper, we propose a coupled porous-free flow and geomechanics model to study the multi-field coupling mechanism in fractured-vuggy reservoirs. Darcy's law is used to describe transport in matrix and fractures, and the Navier–Stokes equation is used to describe free flow in vugs and large fractures. Linear elastic deformation of porous media, as well as deformation of fractures and vugs, are also considered in the model. A finite element numerical simulation method is adopted to solve the double coupled system by MATLAB and COMSOL. The coupling properties, including porous-free flow coupling, and porous flow-geomechanics coupling are investigated by numerical experiments. Moreover, we systematically analyze how the double coupling behaviors are affected by vugs through parameter sensitivity studies in fractured-vuggy carbonate reservoirs. The results have shown that vugs play a dominant role in the production compared with fractures. Fluid flow and pressure drop propagation are greatly affected due to the serious cavern storage effects and stress concentration, as well as significant stress interference between vugs. Vug parameters, including number, size, and distribution pattern, dominate the propagation rate of formation pressure drop and stress distribution.