It can be shown by calculation of divergence of a field like $(\frac{x}{\sqrt{x^2 + y^2}}, \frac{y}{\sqrt{x^2 + y^2}})$ that it's divergence is positive. But I can't understand the geometrical essence of this statement.
It doesn't seem that inward flow is less than outward (as the speed is the same). I know it can be interpreted like tendency to escape neighborhood of a point spreading out in different directions.
Can you please explain me what spreading out has in common with flow through a infinitesimal sphere around the point?
A field of constant length and direction certainly represents "incompressible" flow, but constant length alone does not, roughly because the values of the field may be "spreading out" at each point (as they are for your field $F$), or "squeezing together" at each point (as they are for the field $-F$, which also has constant length). It may help to think of a crowd of people, each walking at the same constant speed, but either leaving or entering a stadium. The density of people per unit area may well change.
There are at least two of ways to get at your question with calculus: one using Gauss's theorem, and one using the definition of a flow.
Let $F = (F_{1}, F_{2})$ be a smooth vector field in the plane, $D$ a small disk centered at a point $(a_{1}, a_{2})$, $S = \partial D$ the boundary circle, and $\mathbf{n}$ the outward unit normal field along $S$. If $dA$ denotes the area element in the plane and $ds$ denotes the arc length element along $S$, Gauss's theorem asserts that $$ \iint_{D} \operatorname{div}(F)\, dA = \int_{\partial D} (F \cdot \mathbf{n})\, ds. $$ Since $F$ is smooth, the left-hand side is approximately equal to $$ \operatorname{div}(F)(a_{1}, a_{2}) \cdot \operatorname{Area}(D), $$ while the right-hand side is the net flux of $F$ across $S = \partial D$. In the limit as $\operatorname{Area}(D) \to 0$, the net flux across $S$ per unit area enclosed is the divergence of $F$ at the center of $D$.
With notation as above, and working locally in space and time, there exists a flow map $\Phi:D \times (-\varepsilon, \varepsilon) \to \mathbf{R}^{2}$ satisfying $$ \frac{\partial \Phi}{\partial t}(x_{1}, x_{2}, t) = F(x_{1}, x_{2}),\quad \Phi(x_{1}, x_{2}, 0) = (x_{1}, x_{2}). $$ The flow map can be expanded as a series in $t$: $$ \Phi(x_{1}, x_{2}, t) = \left[\begin{array}{@{}l@{}} x_{1} + t F_{1}(x_{1}, x_{2}) + O(t^{2}) \\ x_{2} + t F_{2}(x_{1}, x_{2}) + O(t^{2}) \\ \end{array}\right]. $$ The infinitesimal spatial behavior of the flow map is given by the spatial differential $$ D\Phi(x_{1}, x_{2}, t) = \left[\begin{array}{@{}rr@{}} 1 + t \dfrac{\partial F_{1}}{\partial x_{1}}(x_{1}, x_{2}) + O(t^{2}) & t \dfrac{\partial F_{1}}{\partial x_{2}}(x_{1}, x_{2}) + O(t^{2}) \\ t \dfrac{\partial F_{2}}{\partial x_{1}}(x_{1}, x_{2}) + O(t^{2}) & 1 + t \dfrac{\partial F_{2}}{\partial x_{2}}(x_{1}, x_{2}) + O(t^{2}) \\ \end{array}\right]. $$ Particularly, the area of the image of $D$ under the time $t$ flow is approximately the area of $D$ multiplied by $$ \det\bigl(D\Phi(a_{1}, a_{2}, t)\bigr) = 1 + t \operatorname{div}(F)(a_{1}, a_{2}) + O(t^{2}). $$ Again, $\operatorname{div}(F)$ measures the infinitesimal (in space) rate of change of area with respect to time under the flow of $F$.
(Both explanations generalize with obvious modifications to vector fields on $\mathbf{R}^{n}$. The second can be modified to show that the curl of $F$ measures the axis and angular speed of the flow of a vector field on $\mathbf{R}^{3}$.)