"In general, when we have more than a single pole, i.e a double pole, triple pole etc. is it best to just try and find the Laurent expansion about the pole?"
In my experience, for small orders, usually it's easier to use a well known formula.
Given an holomorphic function $\varphi$ with a pole order $k$ at $z_0$, one has, around $z_0$,
$$\varphi (z)=\sum \limits_{n=-k}^\infty \left(c_n(z-z_0)^n\right),$$ for a certain complex sequence $\left(c_n\right)_{n\in \mathbb Z}$, where $\forall n\in \mathbb Z(n<-k\implies c_n=0)$.
Multiply by $(z-z_0)^k$ to get
$$\begin{align}(z-z_0)^k\varphi(z)&=\sum \limits_{n=-k}^\infty \left(c_n(z-z_0)^{n+k}\right)\\ &=c_{-k}+c_{-k+1}(z-z_0)+\ldots +\underbrace{c_{-1}}_{\text{Res}(\varphi , z_0)}(z-z_0)^{-1+k}+\ldots .\end{align}$$
Differentiating $k-1$ times yields $$\dfrac{\mathrm d^{k-1}}{\mathrm dz^{k-1}}\left(z\mapsto (z-z_0)^k\varphi(z)\right)(z)=(k-1)!c_{-1}+k!c_0(z-z_0)+\ldots.$$
Thus, evaluating the above at $z_0$ gives you the residue since everything on the RHS except for $(k-1)!c_{-1}$ goes astray. However $\varphi$ is not defined at $z_0$, so you can't really evaluate at $z_0$, that's OK, take the limit at $z_0$ which exists by definition of pole of order $k$.
This justifies the formula $$\text{Res}(\varphi, z_0)=\dfrac 1{(k-1)!}\lim \limits_{z\to z_0}\left[\dfrac{\mathrm d^{k-1}}{\mathrm dz^{k-1}}\left(z\mapsto (z-z_0)^k\varphi(z)\right)(z)\right].$$
For this particular example, I think it's better to compute $c_{-1}$ by using the series expansion.
Set $g\colon \mathbb C\to \mathbb C, z\mapsto \cos(z)-1, S=\{2k\pi i\colon k\in \mathbb Z\}, h\colon \mathbb C\setminus S\to \mathbb C, z\mapsto \dfrac1{\left(e^z-1\right)^2}$, and $f\colon \mathbb C\setminus S\to \mathbb C, z\mapsto g(z)h(z)$.
Let $k\in \mathbb Z$ be such that $k\neq0$ and write $s=2k\pi i$.
'Around' $s$ each of the functions above is (holomorphic, hence) analytic.
Therefore, for some sequences $(a_n)_{n\in\mathbb Z}, (b_n)_{n\in \mathbb Z}, (c_n)_{n\in \mathbb Z}$ and for all $z$ in a sufficiently small punctured disc centered at $s$ one has
$$g(z)=\sum \limits_{n\in \mathbb Z}\left(a_n(z-s)^n\right),$$
$$h(z)=\sum \limits_{n\in \mathbb Z}\left(b_n(z-s)^n\right),$$
$$f(z)=\sum \limits_{n\in \mathbb Z}\left(c_n(z-s)^n\right).$$
You've determined that $s$ is a pole of order $2$ of $f$. The same happens to $h$ since $\lim \limits_{z\to s}\left(\left|(z-s)h(z)\right|\right)=\infty$ and $$\lim \limits_{z\to s}\left((z-s)^2h(z)\right)=\lim \limits_{z\to s}\left(\dfrac{2(z-s)}{2(e^z-1)e^z}\right)=\lim \limits_{z\to s}\left(\dfrac{1}{2e^{2z}-e^z}\right)=1\neq 0 \tag{$\huge ☺$}$$
Thus $$\forall n\in \mathbb Z(n<0\implies a_n=0),$$
$$\forall n\in \mathbb Z(n<2\implies b_n=0),$$
$$\forall n\in \mathbb Z(n<2\implies c_n=0).$$
Therefore $$\sum \limits_{n=-2}^\infty\left(c_n(z-s)^n\right)=f(z)=g(z)h(z)=\sum \limits_{n=0}^\infty\left(a_n(z-s)^n\right)\sum \limits_{n=-2}^\infty\left(b_n(z-s)^n\right).$$
Recall that the goal is to find $c_{-1}$.
In particular and more pictorially one has
$$\begin{align}
g(z)h(z)&=\left(a_0+a_1(z-s)+\ldots\right)\left(b_{-2}(z-s)^{-2}+b_{-1}(z-s)^{-1}+\ldots\right)\\
&=\left(a_0b_{-2}(z-s)^{-2}+(a_0b_{-1}+a_1b_{-2})(z-s)^{-1}+\ldots\right).
\end{align}$$
Consequently $c_{-1}=a_0b_{-1}+a_1b_{-2}$.
Because $g$ is analytic, $a_0$ and $a_1$ are just the coefficients of the taylor series of $g$ around $s$, that is, $$a_0=f(s)=\dfrac{e^{i\cdot 2k\pi i}+e^{-2k\pi i\cdot i}}{2}-1=\cosh(2k\pi)-1$$ and $$a_1=f'(s)=-\sin(s)=-\dfrac{e^{2k\pi i\cdot i}-e^{-2k\pi i\cdot i}}{2i}=\dfrac{e^{2k\pi}-e^{-2k\pi}}{2i}=-i\sinh(2k\pi).$$
Where as $\displaystyle b_{-2}=\dfrac{1}{2\pi i}\int _\gamma \dfrac{h(z)}{(z-s)^{-1}}\mathrm dz$ and $\displaystyle b_{-1}=\dfrac{1}{2\pi i}\int _\gamma h(z)\,\mathrm dz$, where $\gamma$ is an appropriate positively parametrized circle, centered at $s$, going around exactly once.
The integrals can be computed using the residue theorem.
By choosing $\gamma$ with a small enough radius, the only singularity inside the $\text{im}(\gamma)$ is $s$.
One gets $$\int _\gamma \dfrac{z-s}{(e^z-1)^2}\mathrm dz=2\pi i\text{Res}\left(z\mapsto \dfrac{z-s}{(e^z-1)^2}, s\right).$$
The formula given in the first part of this answer, together with $\huge ☺$, yields $\text{Res}\left(z\mapsto \dfrac{z-s}{(e^z-1)^2}, s\right)=1$, thus $b_{-2}=1$.
For the other, since $s$ is a pole of order $2$, one gets
$$\begin{align}
\int _\gamma h(z)\,\mathrm dz=\int _\gamma \dfrac{1}{(e^z-1)^2}\mathrm dz&=2\pi i\text{Res}(h,s)\\
&=2\pi i\lim \limits_{z\to s}\left[\dfrac{\mathrm d}{\mathrm dz}\left(z\mapsto (z-s)^2h(z)\right)(z)\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{\mathrm d}{\mathrm dz}\left(w\mapsto \dfrac{w^2}{(e^{w+2k\pi i}-1)^2}\right)(z)\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{\mathrm d}{\mathrm dz}\left(w\mapsto \dfrac{w^2}{(e^{w}-1)^2}\right)(w)\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w(e^w-1)^2-2w^2(e^w-1)e^w}{(e^w-1)^4}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w(e^{w}-1)-2w^2e^w}{(e^{w}-1)^3}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w\left(w+\frac{w^2}2+\ldots\right)-2w^2(1+w+\ldots)}{(w+\ldots)^3}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w^2+w^3-2w^2-2w^3+\ldots}{w^3+\ldots}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{-w^3+\ldots}{w^3+\ldots}\right]\\
&=-2\pi i,
\end{align}$$
therefore $b_{-1}=-1$.
Thus $\text{Res}(f,s)=c_{-1}=1-\cosh(2k\pi)-i\sinh(2k\pi)$. This agrees with wolfram alpha.
Best Answer
Yes, you are correct, $z=1$ is the only pole and it is of order $2$. As regards the residue evaluation note that $$\frac{2z+1}{(z-1)^2}=\frac{2z-2+3}{(z-1)^2}=\frac{3}{(z-1)^2}+\frac{2}{z-1}.$$ Now recall the definition of residue. Can you take it from here?