  We perform numerical simulations of the Kelvin-Helmholtz instability in the
mid-plane of a protoplanetary disk. A two-dimensional corotating slice in the
azimuthal--vertical plane of the disk is considered where we include the
Coriolis force and the radial advection of the Keplerian rotation flow. Dust
grains, treated as individual particles, move under the influence of friction
with the gas, while the gas is treated as a compressible fluid. The friction
force from the dust grains on the gas leads to a vertical shear in the gas
rotation velocity. As the particles settle around the mid-plane due to gravity,
the shear increases, and eventually the flow becomes unstable to the
Kelvin-Helmholtz instability. The Kelvin-Helmholtz turbulence saturates when
the vertical settling of the dust is balanced by the turbulent diffusion away
from the mid-plane. The azimuthally averaged state of the self-sustained
Kelvin-Helmholtz turbulence is found to have a constant Richardson number in
the region around the mid-plane where the dust-to-gas ratio is significant.
Nevertheless the dust density has a strong non-axisymmetric component. We
identify a powerful clumping mechanism, caused by the dependence of the
rotation velocity of the dust grains on the dust-to-gas ratio, as the source of
the non-axisymmetry. Our simulations confirm recent findings that the critical
Richardson number for Kelvin-Helmholtz instability is around unity or larger,
rather than the classical value of 1/4
