Is there such a thing as band-limited non-linear distortion?
Applying a non-linear function will always introduce harmonics, and mixing non-linear functions with sampled versions of continuous signals does add the wrinkle you note above (where high-frequency harmonics are aliased to low frequencies.)

I can think of a few ways to proceed:


You can use an oversampling factor high enough to capture the extra harmonics (up to some arbitrary precision, e.g. your noise floor),
You can use a "softer" clipping function (see, for example, here) that has harmonics that die out sooner than the hard clipper. This is easier to model, but introduces its own distortion at low frequencies.
Building on the approach you suggested above, interpolate your sampled signal (e.g. using a Lagrange or Chebyshev interpolator) to construct a continuous-time model. Then, apply the hard clipper and low-pass in a simulated continuous-time domain. Sample the result.


You can combine (1) and (2). The third approach is complex, but gives you the best control over how much distortion to admit, and will probably scale better to very high fidelity requirements.

For non-linear functions that admit a series expansion (e.g. Taylor/Maclaurin), you can get a decent intuition for how fast the harmonics decay. The Maclaurin expansion of a function $f(x)$ is:

$$f(x)=\sum_{n=0}^\infty \left[ \frac{f^{(n)}(0)}{n!} x^n\right]$$

In your case, $f(x)$ is the clipping function. (You can't do this with a hard clipper, at least not naively!) If you consider the substitution $x=g(t)$, where $g(t)$ is your input signal, the $x^n$ becomes $g(t)^n$, which you can consider to be the convolution of your input signal with itself $n$ times. Thus, for low-pass signals, the $n$th term of the infinite summation has a bandwidth $n$ times that of your signal. To complete the picture, you need to figure out the amplitude associated with each term and decide how many terms in the summation are relevant.

(With a little thought, you might also be able to use this form directly to approximate the filtered non-linearity. That would require a good series representation for the clipper.)