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$.
Here is a solution more related to Evan's approach in the mentioned textbook. Let $u = u(x,t)$ be a solution of the PDE. Then using integration by parts and that $u_t = \Delta u$ we get
$$\begin{align}\frac{d}{dt} \left(\frac{1}{2} \|u\|_{L^2(U)}^2\right) &= \int_U u_tu\ dx = \int_U u \Delta u\ dx\\&= -\int_U |Du|^2 dx \overset{(\ast)}\leq -
\lambda_1 \|u\|^2_{L^2(U)}\end{align}$$
where $(\ast)$ comes from Rayleigh's Formula
$$\lambda_1 = \underset{\substack {u \in H_0^1 (U)\\ u\neq 0}}\min \frac{B[u,u]}{\|u\|^2_{L^2(U)}} = \underset{\substack{u \in H_0^1 (U)\\ u\neq 0}} \min \frac{\int_U |Du|^2 dx}{\|u\|^2_{L^2(U)}} $$
Now let $\eta (s) = \|u(\cdot,s)\|^2_{L^2(U)}$. Then
$$\frac{d}{ds} \left(\eta(s) e^{2\lambda_1 s}\right) = e^{2\lambda_1 s} (\eta'(s) +2\lambda_1 \eta(s)) \leq 0$$
Integrating from $0$ to $t$ w.r.t. $s$ we obtain
$$\eta(t)e^{2\lambda_1 t} \leq \eta (0)$$
Since $\eta (0) = \|u (\cdot, 0)\|^2_{L^2(U)} = \|g\|^2_{L^2(U)}$ the result follows.
Best Answer
A partial answer.
In the book by Gilbarg and Trudinger it is proved the following (Theorem 8.12, page 186 of the last edition). I simplify the statement a little bit.
$f \in L^2(U)$ and $\varphi \in H^2(U)$ are given. If $L$ is strictly elliptic with smooth coefficients and the boundary of the domain $U$ is regular, then the solution of $$ \begin{aligned} &Lu = f \text{ in } U\\ &u-\varphi \in H^2(U) \end{aligned} $$ satisfies $$ \|u\|_{H^2(U)}^2 \leq C(\|u\|_{L^2(U)}^2+\|f\|_{L^2(U)}^2+\|\varphi\|_{H^2(U)}^2). $$ $C$ depends on $L,U$. The theorem is written without squared norms, but you can switch from one to the other easily.
Using this theorem it can be proved that (this is Evans exercise in the case $L = -\Delta$) if $u \in H^2(U) \cap H_0^1(U)$, then $$ \|u\|_{H^2(U)}^2 \leq C(\|u\|_{L^2(U)}^2+\|-\Delta u\|_{L^2(U)}^2). $$ Just set $f = \Delta u \in L^2(U)$. A solution of $$ \begin{aligned} &\Delta w = f \text{ in } U\\ &w = 0 \text { on } \partial U \end{aligned} $$ must satisfy $$ \|w\|_{H^2(U)}^2 \leq C(\|w\|_{L^2(U)}^2+\|-\Delta w\|_{L^2(U)}^2). $$ but $u$ itself is a solution of that problem. Actually the book by Evans is also a reference for this theorem.
To prove the exercise in general one should prove that $$ \theta\|-\Delta w\|_{L^2(U)}^2 \leq (Lu,-\Delta u), $$ where $\theta$ is $L$ ellipticity constant. This can be easily seen in dimension 1, but I did not do the computation in higher dimension. This last point should be checked.