I think the analytical formula in terms of rho is going to be (1.57/2*pi) * sqrt(E / rho * L^2), since the Roark formula is (1.57/2*pi) * sqrt(A* E * g / w * L^2) and the weight per unit length is w = m * g / L = rho * A * g. With your values for L, A, E, and rho, this gives f1 = 16 kHz. Since A does not appear in the analytical formula, this also explains why changing the area in the Comsol model doesn't change the frequency.