What you need is the following:
Let $v \in C^\infty_c(U)$ and $w\in C^\infty(\bar{U})$, we have
$$ \left(\int_U Dv \cdot Dw ~\mathrm{d}x\right)^2 \leq C \int_U |v|^2 \mathrm{d}x \int_U |D^2 w|^{2}\mathrm{d}x \tag{*}$$
This follows by directly integrating by parts (the boundary terms vanish as $v$ has compact support).
Now, given $u \in H^1_0(U) \cap H^2(U)$, let $v_i \to u$ in $H^1_0$ and $w_i \to u$ in $H^2(U)$ where $v_i \in C^\infty_c(U)$ and $w_i \in C^\infty(\bar{U})$.
By the strong convergence in $H^1_0$ and $H^2$ respectively, we have that for any function $f\in L^2$ we have
$$ \lim_{\ell \to \infty}\int_U \partial_{x^j} v_\ell f \mathrm{d}x = \lim_{\ell \to \infty}\int_{U} \partial_{x_j} (v_\ell - u + u) f \mathrm{d}x = \int_{U} \partial_{x^j} u f \mathrm{d}x + \lim_{\ell\to\infty}\int_{U} (\partial_{x_j}v_\ell - \partial_{x_j}u) f \mathrm{d}x $$
The second term on the RHS tends to zero using Cauchy-Schwarz and the assumed convergence of $v_\ell\to u$. Similarly we also have
$$ \lim_{\ell \to \infty}\int_U \partial_{x_j} w_\ell f \mathrm{d}x = \int_{U} \partial_{x_j} u f \mathrm{d}x $$
So we have that
$$ \int_U |Du|^2 \mathrm{d}x = \lim_{i,j\to \infty} \int_U Dv_i \cdot D w_j ~\mathrm{d}x \leq \lim_{i,j\to\infty} C \|v_i\|_{L^2} \|D^2 w_j\|_{L^2} $$
by (*). Since $v_i \to u$ in $H^1_0$, we also have $v_i \to u$ in $L^2$. Similarly as $w_j \to u$ in $H^2$ we have $D^2w \to D^2 u$ in $L^2$. So the RHS is
$$ \lim_{i,j\to\infty} C \|v_i\|_{L^2} \|D^2 w_j\|_{L^2} = C\|u\|_{L^2} \|D^2 u\|_{L^2}$$
and we have the desired result.
We have that $$\tag{1}\int\nabla u_\epsilon\nabla\varphi+\int\beta(u_\epsilon)\varphi=f\varphi,\ \forall\ \varphi\in H_0^1$$
By taking $\varphi=u_\epsilon$, we conclude $$\tag{2}\int|\nabla u_\epsilon|^2\leq\int fu_\epsilon$$
which implies by Holder inequality that $\|\nabla u_\epsilon\|_2$ is bounded, therefore, we can suppose that $u_\epsilon \rightharpoonup u$ in $H_0^1$. Moreover, we can suppose that $u_\epsilon \to u$ in $L^2$. We combine these convergences with $(2)$ to get $$\tag{3}\int |\nabla u|^2\leq \int fu$$
On the other hand, define $\mathcal{K}=\{w\in H_0^1:\ w\geq 0\}$. As you can easily verify, $\int\beta(u_\epsilon)w\leq 0$ in $\mathcal{K}$, hence, from $(1)$, we conclude that $$\tag{4}\int\nabla u_\epsilon\nabla\varphi\geq \int f\varphi,\ \forall\ \varphi\in \mathcal{K}$$
Again, we use the convergences we have, to conclude from $(4)$ $$\tag{5}\int\nabla u\nabla\varphi\geq\int f\varphi,\ \forall\ \varphi\in\mathcal{K}$$
We conclude from $(3)$ and $(5)$ the desired inequality. To conclude that $u\in\mathcal{K}$ and $u$ is unique, you can argue like this:
I - $u$ is the minimizer of some functional $I:\mathcal{K}\to\mathbb{R}$, what is $I$?
II - In this particular setting, for $u$ to minimizes $I$ is equivalently for $u$ to satisfies the variational inequality.
III - The solution of the optimization problems is unique.
As @RayYang suggested this book is a good one to understand it better, however I would like to point out that the main argument here can be found in any good book of convex analysis.
It is worth to note that what we are proving here is that $-I'(u)\in \mathcal{N}_{\mathcal{K}}(u)$, where $\mathcal{N}_{\mathcal{K}}(u)$ is the normal cone of $u$ with respect to $\mathcal{K}$ and when $I$ is a convex differentialbe function, II is valid, i.e. $-I'(u)\in \mathcal{N}_{\mathcal{K}}(u)$ if and only if $u$ minimizes $I$.
Best Answer
The reason why you need to find a space $W$ satisfying $V\subset\subset W \subset\subset U$ is that you need to smoothly extend the function from set $V$, which takes value 1 to the set close to the boundary of $U$, which takes value $0$. So consider the classical solution of Laplace equation with scale and shift terms included to match the structure.