By a well-known theorem of Whitney, any compact subset of $R^n$ coincides with the level set of a $C^\infty$ function:
Let $f$ be such a function for the boundary of the domain. By Sard's theorem $f$ has a regular value $x$ arbitrarily close to $0$. Then a component of $f^{-1}(x)$ yields the desired approximation. To get a subdomain, make sure that $x$ coincides with a value of $f$ inside the domain.
There is an even easier way of doing this, cover the domain by balls of radius $\epsilon$, or take the union of all balls of radius $\epsilon$ contained in the domain. For each ball let $f\colon R^n\to R$ be the function whose zero set coincides with the boundary of the ball (i.e., the distance function from the center of the ball minus the radius). Now multiply all these functions to obtain a function $F$, and take a level set of $F$ close to $0$.
The second method yields a subdomain which is not only smooth but is algebraic. In both cases here, the approximation is with regard to the Hausdorff distance.