In rotorcrafts, the ground resonance may arise from interaction among lead-lag motion of the rotor blades, fuselage inertia, and stiffness and damping characteristics of the shock absorbers in landing gears. Approximate prediction of the ground resonance is often conducted by linear modeling and analysis of characteristic roots of the system, although the effect of nonlinear quadratic damping was discovered recently by applying harmonic balance method. Double-stage oleo-pneumatic shock absorbers, employed for softer landing and more comfortable taxiing, show very strong nonlinear stiffness characteristics depending on dynamic amplitude, which is not known a priori. In this study, how to deal with this nonlinear pneumatic stiffness problem by the harmonic balance method in an iterative manner is presented for an eight-degree of freedom rotorcraft system and simulation results are discussed in two aspects. One is the necessity of the nonlinear treatment of the dual-stage pneumatic stiffness relative to the linear stiffness approximation. The other is how to choose the ratio of pneumatic pressure in two chambers or location of the breakover point relative to the static equilibrium point.