$\sum_{k =1, k \neq j}^{N-1} \csc^2\left(\pi \frac{k}{N} \right)\csc^2\left(\pi \frac{j-k}{N} \right)=?$

Start from the well known formula \begin{equation} 2^{N-1} \prod _{k=1}^N \left[\cos (x-y)-\cos \left(x+y+\frac{2\pi k}{N}\right)\right]=\cos N(x-y)-\cos N(x+y), \end{equation} take logarithmic derivative \begin{equation} \sum _{k=1}^N \frac{\sin(x-y)}{\cos (x-y)-\cos \left(x+y+\frac{2\pi k}{N}\right)}=\frac{N\sin N(x-y)}{\cos N (x-y)-\cos N(x+y)}, \end{equation} rewrite it as \begin{equation} N+\sum _{k=1}^N \frac{\cos(x-y)+\cos \left(x+y+\frac{2\pi k}{N}\right)}{\cos (x-y)-\cos \left(x+y+\frac{2\pi k}{N}\right)}=\frac{2N\sin N(x-y)}{\cos N (x-y)-\cos N(x+y)}\cot(x-y), \end{equation} and simplify to get \begin{equation} \sum_{k=1}^N\cot\left(x-\frac{\pi k}{N}\right)\cot\left(y-\frac{\pi k}{N}\right)=\frac{2N\cot(x-y)\sin N(x-y)}{\cos N(x-y)-\cos N(x+y)}-N. \end{equation} Then take derivatives wrt $x$ and $y$ $$ \sum_{k=1}^N\csc^2\left(x-\frac{\pi k}{N}\right)\csc^2\left(y-\frac{\pi k}{N}\right)=\frac{N}{\sin ^2(x-y)} \left(\frac{N}{\sin ^2Nx}+\frac{N}{\sin ^2Ny}-\frac{2\sin N(x-y)}{\sin Nx\sin Ny}\,\cot (x-y) \right). $$

Yes, it may be simplified. The text below is not probably the shortest way, but it explains how to calculate many similar sums.

We start with $$\sin^2 x=-\frac14(e^{ix}-e^{-ix})^2=-\frac14e^{-2ix}(e^{2ix}-1)^2.$$ So, denoting $\omega_k=e^{2\pi i k/N}$ for $k=0,\ldots,N-1$ we get $$S:=\sum_{1\leqslant k\leqslant N-1,k\ne j}\csc^2\left(\pi \frac{k}{N} \right)\csc^2\left(\pi \frac{j-k}{N} \right)= 16\sum_{k\ne 0,j}\frac{\omega_k \omega_{k-j}}{(\omega_k-1)^2(\omega_{k-j}-1)^2}=\\ 16\omega_j \sum_{k\ne 0,j}\frac{\omega_k^2}{(\omega_k-1)^2(\omega_{k}-\omega_j)^2}.$$ Consider the rational function $f(x):=\frac{x^2}{(x-1)^2(x-a)^2}$ (where $a=\omega_j$). We have $$f(x)= \frac{a^2}{(a - 1)^2 (x - a)^2} + \frac{2 a}{(a - 1)^3 (x - 1)} - \frac{2 a}{(a - 1)^3 (x - a)} + \frac1{(a - 1)^2 (x - 1)^2}.$$ Each of four summands may be easily summed up over $x\in \{\omega_0,\ldots,\omega_{N-1}\}\setminus \{1,a\}$. Indeed, we have $$ \sum_{k} \frac 1{t-\omega_k}=\frac{(t^n-1)'}{(t^n-1)}=\frac{nt^{n-1}}{t^n-1},\,\,\, \sum_{k} \frac 1{(t-\omega_k)^2}= -\left(\frac{nt^{n-1}}{t^n-1}\right)'=n\frac{t^{n-2}(t^n+n-1)}{(t^n-1)^2}. $$ If we want to substitute $t=\omega_j$, we get something like $$ \sum_{k:k\ne j} \frac1{\omega_j-\omega_k}=\left(\frac{nt^{n-1}}{t^n-1}-\frac1{t-\omega_j}\right)_{t=\omega_j}=\left(\frac{(n-1)t^{n}-n\omega_j t^{n-1}+1}{(t^n-1)(t-\omega_j)}\right)_{t=\omega_j}=(n-1)\omega_j^{-1} $$ by L'Hôpital, analogously for the squares.