A large eddy simulation (LES) approach is introduced to enable the study of the nonlinear growth of ballooning modes in a heliotron-type device, by solving fully 3D two-fluid magnetohydrodynamic (MHD) equations numerically over a wide range of parameter space, keeping computational costs as low as possible. A model to substitute the influence of scales smaller than the grid size, at sub-grid scale (SGS), and at the scales larger than it—grid scale (GS)—has been developed for LES. The LESs of two-fluid MHD equations with SGS models have successfully reproduced the growth of the ballooning modes in the GS and nonlinear saturation. The numerical results show the importance of SGS effects on the GS components, or the effects of turbulent fluctuation at small scales in low-wavenumber unstable modes, over the course of the nonlinear saturation process. The results also show the usefulness of the LES approach in studying instability in a heliotron device. It is shown through a parameter survey over many SGS model coefficients that turbulent small-scale components in experiments can contribute to keeping the plasma core pressure from totally collapsing.