A stabilized incompressible smoothed particle hydrodynamics (ISPH) method with the addition of a density invariant relaxation condition in the pressure calculations is applied to simulations of highly nonlinear liquid sloshing problems. By applying the Neumann boundary condition when solving pressure, the performance of the present ISPH method is enhanced significantly. Two large-amplitude free sloshing problems under a resonance sway excitation were carried out in a square and a rectangular tank with filling-depths ratios of 20% and 50% of tank height, respectively, and compared with the available published experimental results. To extend the validation of the method, numerical simulations for sloshing problems with the varying density of a floating body as well as a middle baffle, which also generates strongly nonlinear free surface flow, were conducted. The results showed that the present ISPH method produces smooth pressure distribution and significantly reduces spurious oscillation. The proposed ISPH method was shown to be robust and accurate in long time simulation of highly nonlinear sloshing problems.