Is this the correct solution for finding the samples of phase delay at a given frequency of an IIR BiQuad filter?
You almost got it, but the final result is wrong because the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator, respectively, not the ratio of the phases:

$$z=\frac{u}{v}=\frac{|u|e^{j\phi_u}}{|v|e^{j\phi_v}}=\frac{|u|}{|v|}e^{j(\phi_u-\phi_v)}\Longrightarrow \phi_z=\phi_u-\phi_v$$

Furthermore, for finding the phase delay, you need to divide the phase by the (normalized) frequency in radians (as shown in this answer). The result is the phase delay in samples.