For a rigorous treatment you will want to rely on Alt-Caffarelli, but the formal steps that produce the Euler-Lagrange equations are as follows:

Introduce an infinitesimal variation $u(x)\mapsto u(x)+\epsilon(x)$, and compute the variation of the functional to first order in $\epsilon$. Denote by $V$ the region in $\Omega$ within which $u\geq 0$, and $S$ the surface boundary of $V$ on which $u=0$. The variation is
$$L[u+\epsilon]=\int_V|\nabla (u+\epsilon)|^2\,dx + \int_\Omega\chi_{\{u+\epsilon > 0\}}\,dx$$
$$=L[u]+\int_V 2(\nabla u)\cdot(\nabla\epsilon) \,dx - \int_{S} \frac{\epsilon}{n\cdot\nabla u}\,ds+{\cal O}(\epsilon^2)$$
$$=L[u]-2\int_V \epsilon\nabla^2 u \,dx +\int_{S} \epsilon\left(2n\cdot\nabla  u-\frac{1}{n\cdot\nabla u}\right)\,ds+{\cal O}(\epsilon^2).$$
On the second line I used that $n\cdot\nabla u<0$ for $n$ a unit vector normal to the surface $S$ and pointing outward. On the third line I carry out a partial integration of the volume integral.    
The variation of $L$ vanishes for arbitrary $\epsilon$ if the integrand of the volume integral $\int_V$ and the integrand of the surface integral $\int_{S}$ each vanish identically, so 
$$\nabla^2 u=0\;\;\text{in}\;\;V,$$
$$2(n\cdot\nabla u)^2-1=0\;\;\text{on}\;\;S.$$

Since $\nabla u$ has only a normal component on $S$ (because $u=0$ on $S$), we may equivalently write the boundary term as $2|\nabla u|^2=1$ on $S$.

---

More generally, following the same steps to minimize
$$L[u]=\int_V\left( f(x)|\nabla u|^2+g(x)\chi_{u>0}\right)\,dx$$
gives the variation
$$\delta L[u]=-2\int_V\epsilon\nabla\cdot(f(x)\nabla u)\,dx+\int_S\epsilon\left(2f(x)(n\cdot\nabla u)-\frac{g(u)}{n\cdot\nabla u}\right)\,ds$$
hence the Euler-Lagrange equations
$$\nabla\cdot(f(x)\nabla u)=0\;\;\text{in}\;\;V,$$
$$2f(x)(n\cdot\nabla u)^2=g(x)\;\;\text{on}\;\;S.$$