In the field of cardiac modeling, the mechanical action of the heart is often simulated using finite element methods. These simulations are becoming increasingly challenging as the computational domain is customized to a patient's anatomy, within which large heterogeneous tension gradients are generated via biophysical cell models which drive simulations of the cardiac pump cycle. The convergence of nonlinear solvers in simulations of large deformation mechanics depends on many factors. When extreme stress or irregular deformations are modeled, commonly used numerical methods can often fail to find a solution, which can prevent investigation of interesting parameter variations or use of models in a clinical context with high standards for robustness. This paper outlines a novel numerical method that is straightforward to implement and which significantly improves the stability of these simulations. The method involves adding a compressibility penalty to the standard incompressible formulation of large deformation mechanics. We compare the method's performance when used with both a direct discretization of the equations for incompressible solid mechanics, as well as the formulation based on an isochoric/deviatoric split of the deformation gradient. The addition of this penalty decreases the tendency for solutions to deviate from the incompressibility constraint, and significantly improves the ability of the Newton solver to find a solution. Additionally, our method maintains the expected order of convergence under mesh refinement, has nearly identical solutions for the pressure-volume relations, and stabilizes the solver to allow challenging simulations of both diastolic and systolic function on personalized patient geometries.