Firstly, turning to the $\mathbb{R}^1$ case, we note that $$ -u_{xx}+u^2u_x=-\partial_x\left(u_x-\frac{1}{3}u^3\right)=f(x) $$ and so we get a Chini equation (see Inhomogeneous Bernoulli Equation): $$ -u_x+\frac{1}{3}u^3=\int_{x_0}^xf(x')dx'+C $$ and a general solution to this equation is not known. Turning to $\mathbb{R}^3$, you can see that your equation can be written down in the form $$ \nabla\cdot {\bf v}=-f $$ being $$ {\bf v}=\left(u_x-\frac{1}{3}u^3,u_y-\frac{1}{3}u^3,u_z-\frac{1}{3}u^3\right) $$ and you have to satisfy the Gauss identity $$ \int_{\partial\Omega}{\bf v}\cdot d{\bf S}=-\int_{\Omega}dV f. $$ This means that one has to find three functions $v_1,\ v_2,\ v_3$ that satisfy three different Chini equations. Again, a general solution misses in this case. The conclusion to draw is that you need to specialize the problem to see if this becomes manageable in some case or resort to some perturbation technique.

Firstly, turning to the $\mathbb{R}^1$ case, we note that $$ -u_{xx}+u^2u_x=-\partial_x\left(u_x-\frac{1}{3}u^3\right)=f(x) $$ and so we get a Chini equation (see Inhomogeneous Bernoulli Equation): $$ -u_x+\frac{1}{3}u^3=\int_{x_0}^xf(x')dx'+C $$ and a general solution to this equation is not known. Turning to $\mathbb{R}^3$, you can see that your equation can be written down in the form $$ \nabla\cdot {\bf v}=-f $$ being $$ {\bf v}=\left(u_x-\frac{1}{3}u^3,u_y-\frac{1}{3}u^3,u_z-\frac{1}{3}u^3\right) $$ and you have to satisfy the Gauss identity $$ \int_{\partial\Omega}{\bf v}\cdot d{\bf S}=-\int_{\Omega}dV f. $$ This means that one has to find three functions $v_1,\ v_2,\ v_3$ that satisfy three different Chini equations. Again, a general solution misses in this case. The conclusion to draw is that you need to specialize the problem to see if this becomes manageable in some case or resort to some perturbation technique.