In my research I encountered the following nonlinear PDE$$\nabla^2 f + \frac{1}{2} (\nabla f)^2 = 0,$$which is just the Laplace equation with additional term.
Does this equation belong to a class of well-known nonlinear PDEs? Is there a simple way to solve this equation?
I have found a solution to the one dimensional version which may help you. The context of your research would also be helpful.
The one dimensional version is:
$$\frac{d^2 f}{d x^2} + \frac{1}{2}\bigg(\frac{df}{dx}\bigg)^2 = 0$$
A reduction of order can be accomplished by definining $u = \frac{df}{dx}$. The equation then becomes first order nonlinear.
$$\frac{du}{dx} + \frac{1}{2}u^2 = 0$$
$$\frac{2du}{u^2} = -dx$$
$$\frac{2}{u} = x + C$$
$C$ is just an arbitrary integration constant.
Thus
$$u = \frac{1}{\frac{1}{2}(x+C)}$$
From before:
$$\frac{df}{dx} = \frac{1}{\frac{1}{2}(x+C)} $$
$$f = 2 \ln(\frac{1}{2}|x+C|) + D$$
Where $D$ is the second arbitrary constant. You can plug this into the 1D equation and check that it works.
Now if we look at the full equation:
$$\frac{d^2 f}{d x^2} + \frac{d^2 f}{d y^2}+\frac{d^2 f}{d z^2} + \frac{1}{2}\bigg(\bigg(\frac{df}{dx}\bigg)^2 +\bigg(\frac{df}{dy}\bigg)^2 + \bigg(\frac{df}{dz}\bigg)^2 \bigg)= 0$$
Separation of variables does not work in this form, but if we approximate the squares of the derivatives by just the derivatives themselves ($q^2 \approx q$) then it is.
$$\frac{d^2 f}{d x^2} + \frac{d^2 f}{d y^2}+\frac{d^2 f}{d z^2} + \frac{1}{2}\bigg(\bigg(\frac{df}{dx}\bigg) +\bigg(\frac{df}{dy}\bigg) + \bigg(\frac{df}{dz}\bigg) \bigg)\approx 0$$
Then using the anzatz: $f(x,y,z) = X(x)Y(y)Z(z)$ we get three decoupled equations whose sum must be 0.
$$\frac{d^2 f}{d x_i^2} + \frac{1}{2}\bigg(\frac{df}{dx_i}\bigg) = \alpha_i$$
with
$$\sum_i \alpha_i = 0$$
$i$ just runs through $x$, $y$, and $z$.
I don't feel like writing the work right now, plus you may not really need this solution, but using the same reduction of order method as above ($u_x = \frac{dX}{dx}$), one can find that the solution for $X(x)$ is
$$X(x) = -2C_x e^{-x} + 2\alpha_x x + D_x$$
Where $C_x$ and $D_x$ are the arbitrary integration constants.
Thus the full solution would be:
$$f \approx \prod_i (-2C_ie^{-x_i} + 2\alpha_i x_i + D_i)$$
Which is just expressing $f = XYZ$. This could be plotted using some kind of 3D colormap if you wanted to. I don't know how to quickly at the moment. A simple thing to plot (just by using WolframAlpha) is $-e^{-x} + x$, which looks similar to the solution of the 1D problem.
I hope this helps you find your solution.