For the case that $D$ is bounded, as Leo Moos has noted, more details can be found in a textbook of minimal surface, cf. Lemma 1.26 in T. Colding, W. Minicozzi, *A Course in Minimal Surfaces*.

That is, let $F(X)=\frac{X}{\sqrt{1+|X|^2}},$ then
\begin{equation}
F(\nabla u_1)-F(\nabla u_2)=\left(\int_0^1dF\big(\nabla u_2+t(\nabla u_1-\nabla u_2)\big)
dt\right)(\nabla u_1-\nabla u_2).
\end{equation}
From this, one can conclude that $v=u_1-u_2$ satisfies an equation of the form
$$\operatorname{div}(a_{i,j}\nabla v)\leq 0,$$
where the matrix is defined as $(a_{i,j})=\int_0^1dF(\nabla u_2+t(\nabla u_1-\nabla u_2))
dt.$

In particular, for a unit vector $V$ and a vector $X$, we have
$$dF(x)V=\frac{V}{\sqrt{1+|X|^2}}-\frac{\langle V,X\rangle}{(1+|X|^2)^{\frac{3}{2}}}X.$$
Thus,
$$(1+|X|^2)^{\frac{3}{2}}\langle V, dF(X)V\rangle=(1+|X|^2)-\langle V,X\rangle^2\geq 1, $$
which means $(a_{i,j})$ is positive, therefore, the usual comparison principle gives the claim.

linearPDE. The coefficients of this PDE will depend on $u_1,u_2$, but that's no problem. If I am not mistaken the '$c$-term' should be zero, so you can just apply your normal comparison principles. $\endgroup$