$$(\mathcal{P}_{\varepsilon}) : \left\{\begin{array}{ll} \displaystyle -div\left(A(x)\nabla u_\varepsilon(x)\right)= \dfrac{a(x)}{|u_\varepsilon(x)|+\varepsilon} &\mbox{ in }\Omega \\\\ u_\varepsilon = 0 &\mbox{ over }\partial \Omega \end{array} \right.$$ Assume that $\Omega \in C^{1,1}$, be a bounded open set, $a \in L^\infty, a\geq0 \mbox{ and } a\neq 0$ for $x \in \Omega, A(x)=(a_{ij}(x))$ a uniformly coercive matrix with coefficients $a_{i,j} \in C^{0,1}(\bar{\Omega}$). Then, $\mathcal{P}_\varepsilon$ have a unique solution $u_\varepsilon \in W^1_0\cap W^{1,p}(\Omega)$. Moroever, $\exists\ C_\Omega$ uniformly independent of $\varepsilon$ such that \begin{equation} \| \nabla u_\varepsilon \|_{L^2} \leqslant C_\Omega \int_\Omega a(x)\,dx. \end{equation}
I'm not sure how to get this inequality, is there any help?
I am a bit rusty on this, so work with me. We are not proving existence or uniqueness, only answering from where the inequality might come. Since $A(x)$ is uniformly coercive, that means that since
$${\rm -div}(A(x)\nabla u_\epsilon)\ge0,$$ then we can apply maximum principle to $u_\epsilon$ (is super-harmonic). Therefore, the minimum of $u_\epsilon$ is attained in $\partial\Omega$, where it is zero. Therefore $u_\epsilon>0$ in $\Omega$ and $|u_\epsilon|=u_\epsilon$. We can re-write the original problem as $$ -(u_\epsilon+\epsilon)\,{\rm div}(A\nabla u_\epsilon)=a~~.\quad(1)$$ Since $$f\, {\rm div}(A\nabla f)+(\nabla f)^T A (\nabla f)={\rm div}(fA\nabla f)$$ we have the following generalized Green's identity $$\int_\Omega f {\rm div}(A\nabla f)+(\nabla f)^T A (\nabla f)dx=\int_{\partial\Omega}f(A\nabla f)\cdot d\vec{S}.$$ We apply this identity to Eq.(1) and we get, using $u|_{\partial\Omega}=0$, \begin{align} \int_\omega (\nabla u_\epsilon)^T A (\nabla u_\epsilon)-\epsilon{\rm div}(A(x)\nabla u_\epsilon)dx&=\int_\Omega a(x)dx\\ \int_\omega (\nabla u_\epsilon)^T A (\nabla u_\epsilon)&=\int_\Omega a(x)\left(1-\frac{\epsilon}{\epsilon+u_\epsilon}\right)dx\\&\le\int_\Omega a(x)dx~~, \end{align} where in the second step we use that $u_\epsilon$ is solution. But $\alpha|v|^2\le|v^TAv|\le \beta |v|^2$ for two constants $0<\alpha\le\beta$ in general since $A$ components are continuous in $\bar{\Omega}$ and form coercive matrix. Therefore $$ ||\nabla u_\epsilon||^2_2\le\frac{1}{\alpha}\int_\Omega a(x)dx, $$ where constant $\alpha$ is uniform.