In order to keep an almost constant number of neighbors contributing at each particle, use
smoothing length varying in time and in space.
Consider
${d}_{i}$
the smoothing length related to particle
$i$
;
(1)
${W}_{j}\left(i\right)=\widehat{W}\left({x}_{i}{x}_{j},\frac{{d}_{i}+{d}_{j}}{2}\right)$
and
$\nabla {W}_{j}\left(i\right)=grad{}_{{x}_{i}}\left[\widehat{W}\left(x{x}_{j},\frac{{d}_{i}+{d}_{j}}{2}\right)\right]$
if kernel correction
or
(2)
${W}_{j}\left(i\right)=W\left({x}_{i}{x}_{j},\frac{{d}_{i}+{d}_{j}}{2}\right)$
and
$\nabla {W}_{j}\left(i\right)=grad{}_{{x}_{i}}\left[W\left(x{x}_{j},\frac{{d}_{i}+{d}_{j}}{2}\right)\right]$
without kernel correction
At each time step, density is updated for each particle
$i$
, according to:
(3)
$${\frac{d\rho}{dt}}_{i}={\rho}_{i}\nabla \cdot {v}_{i}$$
with
(4)
$$\nabla \cdot {v}_{i}={\displaystyle \sum \frac{{m}_{j}}{{\rho}_{j}}\left({v}_{i}{v}_{j}\right)}\cdot \nabla {W}_{j}\left(i\right)$$
Where,

${m}_{j}$
 Mass of a particle
$i$

${\rho}_{i}$
 Density

${v}_{j}$
 Velocity
Strain tensor is obtained by the same way when non pure hydrodynamic laws are used or in
the other words when law uses deviatoric terms of the strain tensor:
(5)
$${\frac{d{v}^{\alpha}}{d{x}^{\beta}}}_{i}={\displaystyle \sum \frac{{m}_{j}}{{\rho}_{j}}\left({v}_{i}^{\alpha}{v}_{j}^{\alpha}\right)}\frac{d{W}_{j}}{d{x}^{\beta}}\left(i\right),\alpha =\mathrm{1...3},\beta =\mathrm{1...3.}$$
Next the constitutive law is integrated for each particle. Then Forces are computed
according to:
(6)
$${{m}_{i}\frac{dv}{dt}}_{i}={\displaystyle \sum _{j}{V}_{i}}{V}_{j}\left[{p}_{i}\nabla {W}_{j}\left(i\right){p}_{j}\nabla {W}_{i}\left(j\right)\right]{\displaystyle \sum _{j}{m}_{i}}{m}_{j}{\pi}_{ij}\frac{\left[\nabla {W}_{j}\left(i\right)\nabla {W}_{i}\left(j\right)\right]}{2}$$
Where
${\rho}_{i}$
and
${p}_{j}$
are pressures at particles
$i$
and
$j$
, and
${\pi}_{ij}$
is a term for artificial viscosity. The expression is more
complex for non pure hydrodynamic laws.
Note: The previous equation reduces to the following
one when there is no kernel correction:
(7)
$${{m}_{i}\frac{dv}{dt}}_{i}={\displaystyle \sum _{j}{V}_{i}}{V}_{j}\left[{p}_{i}+{p}_{j}\right]\nabla {W}_{j}\left(i\right){\displaystyle \sum _{j}{m}_{i}}{m}_{j}{\pi}_{ij}\nabla {W}_{j}\left(i\right),$$
since
$\nabla {W}_{i}\left(j\right)=\nabla {W}_{j}\left(i\right)$
Then, in order particles to keep almost a constant number of neighbors into their kernels (
$\rho {d}^{3}$
is kept constant), search distances are updated according
to:
(8)
$$\frac{d({d}_{i})}{dt}={{d}_{i}\frac{\nabla \cdot {v}_{i}}{3}}_{}$$