Obtaining the Green's function for a 2D Poisson equation ( in polar coordinates) I am trying to solve the following BVP within an annular region of radii $r_1$, and $r_2$ :
$$
\begin{cases}
\nabla^2u=f\\
u(r_1) = p\\ u(r_2) = q
\end{cases}
$$
If we define an auxiliary problem in terms of Green's function as
$$
\begin{cases}
\nabla^2G=\delta^2(r-r')\\ 
G(r_1) = 0 \\ 
G(r_2) = 0\\
\end{cases}
$$
We have the solution of u (as given by Green's identities as)
$$\DeclareMathOperator{\Dm}{\operatorname{d}\!}
u = \oint u  \frac{\partial G}{\partial n} \Dm S + \int Gf \Dm V   \tag{Eqn. A}
$$
How do I proceed to obtain the form of the Green's function ?
I understand that G for a finite boundary problem is done by superposition :
$G = G_{Freespace} + G_{Homogeneous}$
From my little searching I found that,
$G_{Freespace} = Aln(r-r')$, and
$G_{Homogeneous}=a_0 + a_nr^ncos(n\phi) + b_nr^nsin(n\phi)$
However, the expected solution ( from a paper) that I am seeking is of the form:
$2\pi G=H_0 (r_1,r_2) + \sum_{n=1}^{\infty} H_n(r_1,r_2) cos(n(\phi_1 - \phi_2))$
The form of $H_0$ and $H_n$ are given in the attachment below.
]1
How do I obtain the solution above from the problem posed ?
 A: "The Green's function" is a (very common) grammatical error that Clifford Truesdell I think once called an atrocity to the English language. (Complains a German... ;) ) But you're excused, the paper you cite does it too.^^
Edit: My brain has betrayed me, I thought to have read this in Truesdell's "An iditiot's fugitive essays on Science", but it's actually a footnote in Rohrlich's "Classical Charged Particles", Chap. 4-7. Full quote:

The commonly used expressions "the Green's function" and "a Green's function" represent an atrocity to the English language. I doubt that those who use them ever refer to "a Shakespeare's sonnet".

More on-topic: You're halfway there! Usual technique is to use the free-space and homogeneous contributions you give (the latter is missing a sum!) and try to match the boundary conditions. Here you need to fulfill $G=0$ at $r_1$ and $r_2$. Have fun comparing coefficients and using trigonometric identities. I might have a closer look later myself.
Edit: Some more info to get into the solution process.
The solution of problems like this in cylindrical coordinates is a common problem, and you could have a look at J.D. Jackson: Classical Electrodynamics, Chapter 3.11; Panofsky&Phillips: Classical Electricity and Magnetism, Chapter 4-9f or W.A. Strauss: Partial Differential Equations, Chapter 10.2, for example.
To find the Green function as the sum of the free-space and homogeneous conribution, let's start with the free-space contribution:
It reads
$$
G_f\left(\vec{r},\vec{r}'\right) = - 2 \pi \ln\left(\frac{\left|\vec{r}-\vec{r}'\right|}{r_0}\right)
$$
with an arbitrary constant $r_0$ (if this arbitrariness weirds you out, have a look in Panofsky&Phillips, but it's not that important).
Now to express that in cylindrical coordinates with radius $r$/$r'$ and polar angle $\phi$/$\phi'$ ($r_1$, $r_2$ etc. in the paper, which in your notation are the constant inner and outer radii)
$$
G_f\left(\vec{r},\vec{r}'\right) = - 2 \pi \ln\left(\frac{\sqrt{r^2+\left(r'\right)^2-2rr'\cos\left(\phi-\phi'\right)}}{r_0}\right)
$$
by the law of cosines. But we need to do some work on this to obtain the form in the paper. J.D. Jackson might help you here a bit.
For the homogeneous part we use your ansatz. This stems from a separation ansatz (see the cited references) and contains the sine and cosine contributions (which are eigenfunctions of the second-derivative-operator on a circle) and powers of the radius $r$ with coefficients later to be matched to fullfill the radial part of the Poisson equation and the boundary conditions. Generally we could start off with something like
$$
G_h\left(\vec{r},\vec{r}'\right) = \sum_n r^n\left(a_n \sin\left(n\phi\right) + b_n \cos\left(n\phi\right)\right) \sum_m r'^m\left(a'_m \sin\left(m\phi'\right) + b'_m \cos\left(m\phi'\right)\right)\,,
$$
which looks bad, but can be simplified immediately by using symmetries of the problem.
Invariance under rotation tells us right away that the dependence on $\phi$ and $\phi'$ can only be a dependence on their difference, so we'd expect something like
$$
... \left(a_n \sin\left(n\left(\phi-\phi'\right)\right) + b_n \cos\left(n\left(\phi-\phi'\right)\right)\right)\,,
$$
but then further we can use that $G\left(\vec{r},\vec{r}'\right) = G\left(\vec{r}',\vec{r}\right)$ and we're only left with the cosines.
You see the solution in the paper come along, right? Still, way to go.
