Incompressible fluid definition In a fluid mechanics course I found that an incompressible fluid flow means literally:
$$\rho = \text{constant} \quad \forall \vec r,\, \forall t$$
Where $\vec r = (x, y, z)$
In my understanding, this means literally that the fluid density is uniform? (Am I wrong?)
In the other hand we can find also that an incompressible fluid means:
$$\dfrac{D\rho}{Dt} = 0$$ which does not necessarily mean that the fluid density is uniform.
What's wrong here?
 A: The general definition of an incompressible flow is $\frac{D\rho }{Dt}=0$ : the density of a fluid particle does not change along its path.
For example, if $\overrightarrow{v}=v(x)\overrightarrow{{{e}_{x}}}$ and $\rho =\rho (y)$ : the path lines are horizontal lines and on such a line, the density does not change.
The condition $\rho =cst$ is a particular case ("incompressible fluid" rather than "incompressible flow").
But frequently, one mean $\rho =cst$ when one speak of incompressible fluid !
Sorry for my poor english !
A: In short
$\rho = constant$ is referred to as incompressible fluid. It is a property of the fluid itself: The equation of state allows the density to be assumed approximately constant regardless of the precise flow. Generally speaking this is a property of liquids and is particularly useful in hydrostatics.
In a moving fluid static pressure and kinetic energy may be converted, as e.g. stated by Bernoulli's law. A vanishing Lagrangian derivative $\frac{D \rho}{D t} = 0$ in this context is termed incompressible flow and is a property of the particular flow. It allows for density fluctuations but along any particular streamline a fluid parcel cannot be compressed.
If the Mach number of a flow is small $Ma \lesssim 0.3$ and further the density is not changed by changes in temperature and chemical decomposition the equation of state of a perfect gas leads to only small density fluctuations due to pressure changes. Thus for only small temperature changes it asymptotically approaches the behaviour of an incompressible fluid. This is referred to as low Mach number approximation.
Often people incorrectly do not distinguish between the three as the results are similar: The stress tensor simplifies and for an incompressible fluid the energy equation decouples from the momentum equation as well. Nonetheless air is not an incompressible fluid but may behave like it in the limits of incompressible flows, e.g. for low Mach numbers.

I will try to give a more detailed explanation.
Most books on fluid dynamics introduce the concept of incompressibility in the beginning - often before even deriving the Navier-Stokes equations. I will start the other way around - with the full conservation equations - and explain step by step why each of these simplifications might be useful.
The full conservation equations
The fundamental equations of compressible fluid dynamics describe the conservation of mass \eqref{1}, momentum \eqref{2} and energy \eqref{3} on a continuum level.
$$\frac{\partial \rho}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j )}{\partial x_j }=0 \tag{1}\label{1}$$
$$\frac{\partial (\rho u_i )}{\partial t} + \sum\limits_{j \in \mathcal{D}}\frac{\partial (\rho u_i u_j )}{\partial x_j} = \sum\limits_{j \in \mathcal{D}} \frac{\partial \sigma_{ij}}{\partial x_j } + \rho g_i \tag{2}\label{2}$$
$$\frac{\partial (\rho e)}{\partial t} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho u_j e)}{\partial x_j} = - \sum\limits_{j \in \mathcal{D}} \frac{\partial q_j}{\partial x_j} + \sum\limits_{i, j \in \mathcal{D}} \frac{\partial (\sigma _{ij} u_i)}{\partial x_j} + \sum\limits_{j \in \mathcal{D}} \rho u_j g_j \tag{3}\label{3}$$
The total energy is given by the combination of internal $e_{in}$ and macroscopic energy $e := e_{in} + \sum\limits_{j \in \mathcal{D}} \frac{u_j u_j}{2}$ while the local heat flux $q_i$ according to Fourier's law is assumed to be proportional to the gradient of the temperature
$$q_i = - k \frac{\partial T}{\partial x_i}. \tag{4}\label{4}$$
and the stress tensor for an isotropic Newtonian fluid (derivation here) is given by
$$\sigma_{ij} = - p \delta_{ij} + \underbrace{ 2 \mu S_{ij} - \frac{2}{3} \mu \sum\limits_{k \in \mathcal{D}} S_{kk} \delta_{ij} }_{\tau_{ij}} \tag{5}\label{5}$$
with the rate of strain tensor
$$S_{ij} := \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right). $$
When furthermore supplying an equation of state that connects pressure, density and temperature this system of coupled differential equations is closed.
Simplifications
The general conservation equations are far too complex to extract any information out of them analytically. We will have to introduce major simplifications. The first thing we could try to do is derive a system of equations with constant coefficients. The advantage of constant coefficients is the law of similarity: A system with the same dimensionless coefficient in the differential equations describing its behaviour with the same type of boundary conditions and geometry must have the same solution. This insight allows us to rescale parameters in such a way that the dimensionless parameters of interest, that appear in the corresponding differential equations, are not changed and perform model experiments which allow us to make conclusions about systems of different dimensions.
Non-conservative form
We will start off by recombining the conservation equations to derive their non-conservative form which is convenient as it lets us see which terms are not constant.
Let's apply the chain rule to equation \eqref{2} and then subtract \eqref{1}
$$\frac{\partial u_i}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial u_i}{\partial x_j} = - \frac{1}{\rho} \frac{\partial p}{\partial x_i } + \frac{1}{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial \tau_{ij}}{\partial x_j } + g_i. \tag{6}\label{6}$$
Applying the chain rule to \eqref{3} and subtracting equation \eqref{1} we yield the so called non-conservative form of the conservation equations that still holds for all Newtonian fluids
$$\frac{\partial e}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial e}{\partial x_j} = - \frac{1}{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial q_j}{\partial x_j} + \frac{1}{\rho} \sum\limits_{i, j \in \mathcal{D}} \frac{\partial (\sigma _{ij} u_i)}{\partial x_j} + \sum\limits_{j \in \mathcal{D}} u_j g_j. \tag{7}\label{7}$$
Now we subtract the dot product of \eqref{6} and velocity $u_i$ and apply the chain rule to obtain the mechanical energy
$$\frac{1}{2} \sum\limits_{i \in \mathcal{D}} \frac{\partial (u_i u_i)}{\partial t} + \frac{1}{2} \sum\limits_{i, j \in \mathcal{D}} u_j \frac{\partial (u_i u_i )}{\partial x_j} = - \frac{1}{\rho} \sum\limits_{i \in \mathcal{D}} u_i \frac{\partial p}{ \partial x_i } + \frac{1}{\rho} \sum\limits_{i \in \mathcal{D}} u_i \frac{\partial \tau_{ij}}{\partial x_j } + \sum\limits_{i \in \mathcal{D}} u_i g_i \tag{8}\label{8}$$
Applying again the chain rule to \eqref{7} and subtracting \eqref{8}, the expression
$$\frac{D e_{in}}{D t} = \frac{\partial e_{in}}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial e_{in}}{\partial x_j}= - \frac{1}{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial q_j}{\partial x_j} - \frac{p}{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial u_j}{\partial x_j} + \frac{1}{\rho} \sum\limits_{i, j \in \mathcal{D}} \tau _{ij} \frac{\partial u_j}{\partial x_i} \tag{9}\label{9}$$
for the thermal energy can be found.
To simplify things even further we will introduce the specific enthalpy $h := e_i + \frac{p}{\rho}$
$$\frac{D h}{D t} = \frac{D e_{in}}{D t} + \frac{1}{\rho} \frac{D p}{D t} - \frac{p}{\rho^2} \underbrace{ \frac{D \rho}{D t} }_{\frac{\partial \rho}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial \rho}{\partial x_j}} = \frac{D e_{in}}{D t} + \frac{1}{\rho} \frac{D p}{D t} + \frac{p }{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial u_j}{\partial x_j} \tag{10}\label{10}$$
where I used equation \eqref{1} to simplify the last term. Finally combining \eqref{9} and \eqref{10} we end up with
$$\frac{\partial h}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial h}{\partial x_j} = \frac{1}{\rho} \left( \frac{\partial p}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial p}{\partial x_j} \right) - \frac{1}{\rho} \sum\limits_{j \in \mathcal{D}} \frac{\partial q_j}{\partial x_j} + \frac{1}{\rho} \sum\limits_{i, j \in \mathcal{D}} \tau _{ij} \frac{\partial u_j}{\partial x_i}. \tag{11}\label{11}$$
Non-dimensional conservation equations
Now if we want to have constant coefficient in front of every single term we will have to assume constant transport coefficients $k$ and $\mu$. If we do so, we are able to introduce the following dimensionless quantities into equations \eqref{6}, \eqref{7} and \eqref{11}:
$$x_i^*=\frac{x_i}{L}, \phantom{abc} u_i^*=\frac{u_i}{U}, \phantom{abc} \rho^*=\frac{\rho}{\rho_0}, \phantom{abc} T^*=\frac{\Delta T}{\Delta T_0}, \phantom{abc} g_i^*=\frac{g_i}{g}, \phantom{abc} t^*=\frac{t}{\frac{L}{U}},\phantom{abc} p^*=\frac{p}{\rho_0 U^2} \phantom{ab}$$
and furthermore assuming the simplest fluid, a perfect gas, with $h = c_p T$, this results in the dimensionless conservation equations
$$\frac{\partial \rho^*}{\partial t^*} + \sum\limits_{j \in \mathcal{D}} \frac{\partial (\rho^* u_j^* )}{\partial x_j^* }=0$$
$$\rho^* \frac{\partial u_i^*}{\partial t^*} + \rho^* \sum\limits_{j \in \mathcal{D}} u_j^* \frac{\partial u_i^*}{\partial x_j^*} = - \frac{\partial p^*}{ \partial x_i^* } + \frac{1}{Re} \sum\limits_{j \in \mathcal{D}} \frac{\partial \tau_{ij}^*}{\partial x_j^* } + \frac{1}{Fr^2} g_i^*$$
$$\rho^*  \frac{\partial T^*}{\partial t^*} + \rho^* \sum\limits_{j \in \mathcal{D}} u_j^* \frac{\partial T^*}{\partial x_j^*} = Ec \left( \frac{\partial p^*}{\partial t^*} + \sum\limits_{j \in \mathcal{D}} u_j^* \frac{\partial p^*}{\partial x_j^*} \right) + \frac{1}{Pr Re} \sum\limits_{j \in \mathcal{D}} \frac{\partial}{\partial x_j^*} \left( \frac{\partial T^*}{\partial x_j^*} \right) + \frac{Ec}{Re} \sum\limits_{i, j \in \mathcal{D}} \tau _{ij}^* \frac{\partial u_i^*}{\partial x_j^*} \tag{12}\label{12}$$
with the corresponding dimensionless numbers
$Re := \frac{U L}{\nu} \phantom{abc} \frac{\text{inertial forces}}{\text{viscous forces}} \phantom{abc}$ Reynolds number
$\phantom{abc} Ec := \frac{U^2}{c_P \Delta T_0} \phantom{abc} \frac{\text{heat dissipation potential}}{\text{advective transport}} \phantom{abc}$ Eckert number
$Fr := \frac{U}{\sqrt{g \, L}} \phantom{abc} \frac{\text{flow inertia}}{\text{gravity}} \phantom{abc}$ Froude number
$\phantom{abc} Pr := \frac{\mu c_P}{k}=\frac{\nu}{a} \phantom{abc} \frac{\text{viscious diffusion rate}}{\text{thermal diffusion rate}} \phantom{abc}$ Prandtl number
Let's see what more we can simplify with special assumptions about the only leftover material parameter, the density.
Incompressible flow
If we assume Lagrangian derivative of the density to vanish
$$\frac{D \rho}{D t} = \frac{\partial \rho}{\partial t} + \sum\limits_{j \in \mathcal{D}} u_j \frac{\partial \rho}{\partial x_j} = 0 \tag{13}\label{13}$$
we would be left with a very simple condition: \eqref{1} and \eqref{13} can be combined to
$$\sum\limits_{j \in \mathcal{D}} \frac{\partial u_j}{\partial x_j} = 0,$$
the so called divergence free condition. This means that a fluid parcel is not compressed along its way on a stream line as the last term in \eqref{5}, the volumetric dilatation, vanishes. This is a property of the flow and as a benefit we are left with a simpler stress tensor.
Incompressible fluid
An incompressible material on the other hand is a property of the fluid that must be compatible with the equation of state: 
$$ \sum\limits_{j \in \mathcal{D}} \frac{\partial \rho}{\partial x_j} = 0$$
holds - the density is assumed constant throughout the flow field $\rho = constant$. As can be seen from the equation of state of water, the Tait equation
$$p - p_0 = C \left[ \left( \frac{\rho}{\rho_0} \right)^m - 1 \right], $$
(with $C \approx 0.32 GPa$ and $m \approx 7$) this is clearly the case for water: Large changes of pressure are connected to only small changes in density.
The stress tensor is equal to the one of an incompressible flow but furthermore the continuity equation loses its meaning as there is no equation of state coupling pressure and density anymore $\rho \neq \rho(p)$. The energy equation is decoupled from the equation system: Momentum equation can be solved independently and the temperature be calculated as a simple advection-diffusion equation afterwards from the resulting velocity field. From the conservation equations one is able to recover pressure gradients but its absolute value can only be determined with a trick (Poisson's equation).
Low Mach number approximation
Now as can be seen from \eqref{12} the energy equation also decoupled in the limit of $Ec \to 0$. Furthermore one can make use of the Mach number
$$Ma := \frac{U}{c_s} \phantom{distance} \frac{\text{ordered kinetic energy}}{\text{random kinetic energy}}$$
and the speed of sound of a perfect gas $c_s := \left( \frac{\partial p}{\partial \rho} \right)_{S = const} = \sqrt{\gamma R_m T}$ to rewrite the Eckert numbert to
$$Ec=\frac{U^2}{c_P \Delta T_0} \frac{c_s^2}{c_s^2} = Ma^2 (\gamma - 1) \frac{T_0}{\Delta T_0}.$$
Thus for non-isothermal flows $\Delta T_0 \neq 0$ the energy equation will decouple for Mach numbers approaching zero as the Eckert number will be neglibly small as well. From compressible gas dynamics one can derive the density ratio
$$\frac{\rho_0}{\rho} = \left( 1 + \frac{\gamma - 1}{2} Ma^2 \right)^\frac{1}{\gamma - 1}.$$
As can be seen every flow is somewhat compressible but if the Mach number is less than $0.3$ the density changes due to pressure are less than $5\%$ (for realistic heat capacity ratios $1.1 \lesssim \gamma \lesssim 1.8$) and the fluid can be approximated with good accuracy as incompressible: Due to $Ma = \frac{U}{c_s} \ll 1$ one can assume that $c_s^2=\left( \frac{\partial p}{\partial \rho} \right)_S \rightarrow \infty$ and therefore $\left( \frac{\partial \rho}{\partial p} \right)_S \rightarrow 0$. This means this assumption is an asymptotic approximation of a truly incompressible fluid $\rho \neq \rho(p)$ (which itself has an indefinitely large speed of sound) if we neglect changes due to chemical composition or temperature.
Summary
So the point of all these concepts is either a constant density in space and time (incompressible fluid, low Mach number approximation to incompressible fluid for isothermal flows), simpler stress tensor (compressible flow, incompressible fluid, low Mach number approximation to incompressible fluid) or a decoupling of the energy equation from the momentum equation (incompressible fluid, low Mach number approximation to incompressible fluid for isothermal flows). As their manifestation is pretty similar a lot of books do not distinguish them clearly and refer to the conservation equations of an incompressible fluid as "incompressible Navier-Stokes equations".
A: The definition of incompressible is often unclear and changes depending on which community uses it. So let's look at some common definitions:
Constant density
This means the density is constant everywhere in space and time. So:
$$
\frac{D\rho}{Dt} = \frac{\partial \rho}{\partial t} + \vec{u}\cdot\nabla{\rho} = 0
$$
Because density is constant everywhere in space and time, the temporal derivative is zero, and the spatial gradient is zero. 
Low Mach Number
This shows up when the flow velocity is relatively low and so all pressure changes are hydrodynamic (due to velocity motion) rather than thermodynamic. The effect of this is that $\partial \rho / \partial p = 0$. In other words, the small changes in pressure due to flow velocity changes do not change the density. This has a secondary effect -- the speed of sound in the fluid is $\partial p/\partial \rho = \infty$ in this instance. So there is an infinite speed of sound, which makes the equations elliptic in nature. 
Although we assume density is independent of pressure, it is possible for density to change due to changes in temperature or composition if the flow is chemically reacting. This means:
$$
\frac{D\rho}{Dt} \neq 0
$$
because $\rho$ is a function of temperature and composition. If, however, the flow is not reacting or multi-component, you will also get the same equation as the constant density case:
$$
\frac{D\rho}{Dt} = 0
$$
Therefore, incompressible can mean constant density, or it can mean low Mach number, depending on the community and the application. I prefer to be explicit in the difference because I work in the reacting flow world where it matters. But many in the non-reacting flow communities just use incompressible to mean constant density. 
Example of non-constant density
Since it was asked for an example where the material derivative is zero but density is not constant, here goes:
$$
\frac{D\rho}{Dt} = \frac{\partial \rho}{\partial t} + \vec{u}\cdot\nabla \rho = 0
$$
Rearrange this:
$$
\frac{\partial \rho}{\partial t} = -\vec{u}\cdot\nabla\rho
$$
gives a flow where $\rho \neq \text{const.}$ yet $D\rho/Dt = 0$. It has to be an unsteady flow. 
Is there another example of steady flow? In steady flow, the time derivative is zero, so you have:
$$
\vec{u}\cdot\nabla\rho = 0
$$
If velocity is not zero, $\vec{u} \neq 0$, then we have $\nabla \rho = 0$ and so any moving, steady flow without body forces (gravity) or temperature/composition differences must have constant density. 
If velocity is zero, you can have a gradient in density without any issues. Think of a column of the atmosphere for example -- density is higher at the bottom than the top due to gravity, and there is no velocity. So again, $D\rho/Dt = 0$ but density is not constant everywhere. The challenge here of course is that the continuity equation is not sufficient to describe the situation since it becomes $0 = 0$. You would have to include the momentum equation to incorporate the gravity forces. 
