Bilinear Transform (Tustin's Method) applied to the Derivative
It's a matter of perspective if the bilinear transform "fails so miserably" when trying to approximate a derivative. First of all, it's not true that the approximation is quadratic for small $\omega$, it's linear as it should be:

$$D(e^{j\omega})=\frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{-j\omega}}=\frac{2j}{T}\tan\left(\frac{\omega}{2}\right)\approx \frac{j\omega}{T}\quad\textrm{for}\quad |\omega|\ll 1$$

So for small frequencies the discrete-time system behaves as a (discrete-time) differentiator. For larger frequencies we expect the effect of frequency warping, which is exactly what happens here: the infinite continuous-time frequency range is mapped to the interval $\omega\in[0,\pi]$.

In sum, the bilinear transform behaves as expected, but it is also clear that there are better ways to approximate a derivative in discrete-time than applying the bilinear transform to the continuous-time transfer function $H(s)=s$.