  Motivated by the stability of dust laden vortices, in this paper we study the
terminal velocity approximation equations for a gas coupled to a pressureless
dust fluid and present a numerical solver for the equations embedded in the
FARGO3D hydrodynamics code. We show that for protoplanetary discs it is
possible to use the baricenter velocity in the viscous stress tensor, making it
trivial to simulate viscous dusty protoplanetary discs with this model. We also
show that the terminal velocity model breaks down around shocks, becoming
incompatible with the two fluid model it is derived from. Finally we produce a
set of test cases for numerical schemes and demonstrate the performance of our
code on these tests. Our implementation embedded in FARGO3D using an
unconditionally stable explicit integrator is fast, and exhibits the desired
second order spatial convergence for smooth problems.
