  We present a new algorithm for simulating two-fluid gas and dust mixtures in
Smoothed Particle Hydrodynamics (SPH), systematically addressing a number of
key issues including the generalised SPH density estimate in multi-fluid
systems, the consistent treatment of variable smoothing length terms, finite
particle size, time step stability, thermal coupling terms and the choice of
kernel and smoothing length used in the drag operator. We find that using
double-hump shaped kernels improves the accuracy of the drag interpolation by a
factor of several hundred compared to the use of standard SPH bell-shaped
kernels, at no additional computational expense. In order to benchmark our
algorithm, we have developed a comprehensive suite of standardised, simple test
problems for gas and dust mixtures: dustybox, dustywave, dustyshock, dustysedov
and dustydisc, the first three of which have known analytic solutions. We
present the validation of our algorithm against all of these tests. In doing
so, we show that the spatial resolution criterion \Delta < cs ts is a necessary
condition in all gas+dust codes that becomes critical at high drag (i.e. small
stopping time ts) in order to correctly predict the dynamics. Implicit
timestepping and the implementation of realistic astrophysical drag regimes are
addressed in a companion paper.
