By a well-known theorem of Whitney, any closed subset of $R^n$ coincides with the zero 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 when the domain is bounded: cover the domain by balls of radius $\epsilon$, or take the union of alla finite collection of balls of radius $\epsilon$ contained in the domain and covering a compact subset. 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.