It is a very well-known fact that the Dirichlet energy of unit vector fields is ill-defined in the presence of singularities. For my own reference, below I recapitulate two arguments that prove this result.
To first understand why, geometrically, consider any smooth vector field $X\in\Gamma(TM)$. It can be written as $X = \alpha \hat{X}$ where $\hat{X} = \tfrac{X}{|X|}$ is a unit vector field and $\alpha : M \to \mathbb{R}$ is its length. Using the Levi-Civita connection to differentiate $X$, we find \[ \nabla X = (d\alpha)\hat{X} + \alpha\nabla\hat{X}, \] and so the Dirichlet energy of $X$ is \[\begin{align*}E(X) & = \tfrac12\int_{M}|\nabla X|^2~dA \\ & = \tfrac12\int_{M}\big( |d\alpha|^2 + |\alpha|^2|\nabla\hat{X}|^2\big)~dA \\ & = \tfrac12\int_{M}\big( |d\alpha|^2 + |\alpha|^2|\omega|^2\big)~dA,\end{align*} \] where \(\omega\in\Omega^1(M;\mathbb{R})\) is the rotation 1-form of $\hat{X}$ (defined as $\nabla\hat{X} = \omega J\hat{X}$, which can be defined since $\hat{X}$ is a unit vector field implies that any of its variations must be in the orthogonal direction. Now for a unit vector field we have $\alpha \equiv 1$, and so the Dirichlet energy of a unit vector field $\hat{X}$ is the $L^2$-norm of its rotation 1-form. We can evaluate the energy density $|\omega|^2~dA = |\omega(U)|^2 + |\omega(V)|^2$ using any orthonormal basis of tangent vectors $U,V$, and in particular \[|\omega|^2~dA = |\omega(J\hat{X})|^2. \] Around a singularity at a point $p\in M$ the field $\hat{X}$ winds around the circle $k$ times for $k\neq 0$, and so the rotation speed at radius $r > 0$ from the singularity is comparable to $2\pi k / 2\pi r \sim \tfrac{1}{r}$ as $r\to 0^+$, but then \[ E(\hat{X}) \geq \int_{B(p,\varepsilon)} \tfrac12|\nabla \hat{X}|^2~dA \approx \pi \int_{0}^{\varepsilon}\tfrac{1}{r^2}~r\,dr = \infty\] for any $\varepsilon > 0$. One would need to make this estimate precise to get a complete proof.
There is another argument that I also like: first, note that it suffices to prove the result for simply connected, compact, surfaces with boundary with Dirichlet boundary conditions given, by cutting out a smaller piece of $M$ containing a single singularity. Given some boundary data $X_0 \in\Gamma(TM|_{\partial M})$ we can define the Sobolev space \[W^{1,2}_{X_0}(M;TM) = \big\{X\in W^{1,2}(M;TM)~|~X|_{\partial M} = X_0\}, \] where the restriction to the boundary is understood in the trace sense. This is the energy space of the Dirichlet energy in the sense that $E(X)<\infty$ precisely for vector fields in this space. By the direct method of the calculus of variations, we can prove that there exists a minimizer over $W^{1,2}_{X_0}(M;TM)$ of the Dirichlet energy $E$ given that$W^{1,2}_{X_0}(M;TM) \neq \emptyset$. By standard elliptic regularity results, the minimizer is actually of regularity $C^1$, and so this realizes a homotopy between the boundary data $X_0$ and the constant map. Since the degree is homotopy invariant, this implies that $X_0$ must have trivial degree. Therefore, we have shown that if $X_0$ has non-trivial degree then the set $W^{1,2}_{X_0}(M;TM) = \emptyset$ and therefore the Dirichlet energy is infinite.