I assume that $\Omega$ should be open? If so then Green's second identity/integration by parts imply that for any $v \in C^2(\Omega)$,
$$\int_U v\Delta\phi\,dy = \int_U \phi\Delta v\,dy$$
for $U \subseteq \overline{U} \subseteq \Omega$ with smooth enough boundary (say piecewise $C^1$), where $U$ contains the support of $\phi$. Fubini's theorem implies $u_\varepsilon$ is subharmonic, and since $u_\varepsilon$ is smooth, a well known result implies $\Delta u_\varepsilon \geq 0$. One way you can show this is by using Green's identity with the fundamental solution of the Laplacian to find that
$$u_\varepsilon(x) = \frac{1}{\omega_n r^n}\int_{B_r(x)} u_\varepsilon(y)\,dy - C\int_{B_r(x)} \Delta u_\varepsilon(y)\frac{1}{\|x - y\|^{n-2}}\,dy$$
for some positive constant $C$ (I don't remember what $C$ is, but it's not important for the argument and you can find the correct constant by running through it). This part is a fairly standard argument in PDEs; apply Green's second identity to the integral
$$\int_{B_r(x) \setminus B_\delta(x)} \Delta u_\varepsilon(y)\frac{1}{\|x - y\|^{n-2}}\,dy$$
then do a bunch of computations and send $\delta \to 0$. There are also other proofs of this fact elsewhere on the site.
Taking $r \to 0$ and using the fact $u_\varepsilon$ is subharmonic, you can conclude $\Delta u_\varepsilon \geq 0$ (recall the fundamental solution acts like an approximate identity). Therefore
$$\int_U u_\varepsilon\Delta\phi\,dy = \int_U \phi\Delta u_\epsilon\,dy \geq 0$$
for all $\varepsilon > 0$, and the argument is finished upon taking $\varepsilon \to 0$.
