1) Gravity has to be included because the inverted pendulum effect changes the resonant frequencies. The deflection from gravity is tiny but the change in the dynamics is not. The results are not accurate without it. The z-direction probably is unaffected by gravity, but the tilt modes really feel it.
2) You should try a better meshing. Right now COMSOL is calculating a lot of strain/stress in the steel plates. For our purposes, we can imagine that the steel is infinitely stiff. There are options in COMSOL to change the meshing density in the different materials - as we can see from your previous plots, all the action is in the rubber.
3) I don't think the mesh density directly limits the upper measurement frequency. When you redo the swept-sine using the matlab scripting, use a logarithmic frequency grid like we usually do for the Bode plots. The measurement axis should go from 0.1 - 30 Hz and have ~100 points.
In any case, the whole thing looks promising: we've got real solid models and we're on the merge of being able to duplicate numerically the Dugolini-Vass-Weinstein measurements.