We consider the equation $$ -\Delta u =f\chi_\Omega\quad \mbox{in }\mathbb{R}^3, $$ where $\Omega$ is a bounded and smooth domain, $\chi_\Omega$ is the charasteristic funtion of the set $\Omega$ and $f\in L^2(\Omega)$.
Can we deduce the boundary condition: $\nabla u\cdot \nu =0$ on $\partial \Omega$, where $\nu$ is the outer unit normal of $\Omega$?
I have tried to get this result using the weak formulation of the problem, but I didn't get anything like this.
No, you can not guarantee such a thing. Consider for simplicity $\Omega = B(0,1)$.
Now take $u$ depending only on the distance to the center, for example, $$u = -\frac{r^4}{20}+a ,\quad 0\le r\le 1,$$where $a$ is a contants. You can show that $\Delta u = r^2$ inside the unit ball, i.e. is $L^2$.
Now take its continuation on $r>1$ such that $\Delta u=0$: You obtain a Cauchy problem $$\begin{cases}\frac{1}{r^2}\partial_r\left(r^2\partial_r u\right)=0,\quad r>1,\\u(1) = a -\frac{1}{20},\\ u'(1) = -\frac 15.\end{cases}$$ Clearly, this Cauchy problem has a unique solution, hence your $u$ is well-defined on $\Bbb R^3$. On the other hand, it is obvious that the choice of the constant $a$ allows you to get whatever you want for the value of $u'(1)$ (and therefore, for $\nabla u\cdot \nu\big|_{\partial \Omega}$)
edit Note, however, that we don't make any hypothesis on the space where we want ot find $u$. In particular, in my counterexample $u\notin L^2(\mathbb R^3)$.