Thermal inflation models (which feature two inflationary stages) can display damped primordial curvature power spectra on small scales which generate damped matter fluctuations. For a reasonable choice of parameters, thermal inflation models naturally predict a suppression of the matter power spectrum on galactic and sub-galactic scales, mimicking the effect of warm or interacting dark matter. Matter power spectra in these models are also characterised by an excess of power (with respect to the standard ΛCDM power spectrum) just below the suppression scale. By running a suite of N-body simulations we investigate the non-linear growth of structure in models of thermal inflation. We measure the non-linear matter power spectrum and extract halo statistics, such as the halo mass function, and compare these quantities with those predicted in the standard ΛCDM model and in other models with damped matter fluctuations. We find that the thermal inflation models considered here produce measurable differences in the matter power spectrum from ΛCDM at redshifts z > 5 for wavenumbers k ∈ [2, 64] h Mpc −1 , while the halo mass functions are appreciably different at all redshifts in the halo mass range M halo ∈ [10 9 , 10 12 ] h −1 M resolved by our simulations. The halo mass function at z = 0 for thermal inflation displays an enhancement of around ∼ 20% with respect to ΛCDM and a damping at lower halo masses, with the position of the enhancement depending on the value of the free parameter in the model. The enhancement in the halo mass function (with respect to ΛCDM ) increases with redshift, reaching ∼ 40% at z = 5. We also study the accuracy of the analytical Press-Schechter approach, using different filters to smooth the density field, to predict halo statistics for thermal inflation. We find that the predictions with the smooth-k filter we proposed in a separate paper agree with the simulation results over a wider range of halo masses than is the case with other filters commonly used in the literature.