Series expansion for integral including error function $\DeclareMathOperator{\erfc}{erfc}
\DeclareMathOperator{\Ei}{Ei}
$
What is the series expansion of $f$ for small $q$? 
\begin{align}
    U(q) &=  q e^{q^2}\erfc q\\
    I(q,q') &= \int_0^{2\pi} \frac{d\phi}{2\pi} U(\sqrt{q^2 + q'^2 -2qq'\cos\phi})-U(q')\\
    f(q) &= \int_{0}^{\infty} dq'\,I(q,q')
\end{align}
Even better, perhaps you can integrate this analytically to find $f$ [or $I(q,q')]$? When I numerically integrate, the results are not grid dependent for small $q$: i.e. I get consistent results with double the gridpoints and/or double the grid maximum, so it seems to be well-behaved (and the result looks approximately quadratic in $q$).  Here is my attempt to find the series expansion of $f$ by differentiating under both integrals, then integrating over $\phi$ analytically:
\begin{align}
    D_n(q') &\equiv \left[\frac{\partial ^{n} I(q,q')}{\partial q^n}\right]_{q=0}\\
    D_0(q') &= 0\\
    D_1(q') &= \int_0^{2\pi} \frac{d\phi}{2\pi} \left[ \frac{U'(\sqrt{q^2 + q'^2 -2qq'\cos\phi})}{2\sqrt{q^2 + q'^2 -2qq'\cos\phi}} (2q - 2q'\cos\phi) \right]_{q=0}\\
    &= -U'(q') \int_0^{2\pi} \frac{d\phi}{2\pi}   \cos\phi = 0\\
    D_2(q') &=  \int_0^{2\pi} \frac{d\phi}{2\pi} \left[\sin^2\phi \frac{U'(q')}{q'} + \cos^2\phi \, U''(q')\right] \\
    &=  \frac{U'(q')}{2q'} + \frac{U''(q')}2   \\&=  \left( 2q'^3+4q'+\frac1{2q'} \right)e^{q'^2}\erfc q' -\frac{3+2q'^2}{\sqrt\pi} \\
    D_3(q') &= 0\\
    f(q) &= \frac{q^2}2\int_{0}^{\infty} dq'\,D_2(q')  + O(q^4)
\end{align}
$D_2(q') = \frac12 q'^{-1} + O(1)$ so the last integral diverges logarithmically for small $q'$. Also it seems that $D_{2n}(q') = O(q'^{1-2n})$ so the higher terms get worse. 
I don't know why the Taylor series looks like a sum of divergent terms: perhaps the series expansion of $f$ is in non-integer powers of $q$ rather than a Taylor series?
In case it is helpful here is some Mathematica:
U[q_] := q Exp[q^2] Erfc[q]
Dint[qd_, n_] :=  Simplify[Integrate[(D[U[Sqrt[q^2 + qd^2 - 2 q qd Cos[phi]]] -U[qd], {q, n}]) /. q -> 0, {phi, 0, 2 Pi}]/(2 Pi), qd >= 0]

so that Dint[qd, 2] gives the above expression for $D_2(q')$. Integrate[Dint[qd,2], qd] gives an expression in terms of $_{2}F_2$ and that expression is divergent for $q'\to0$. 
EDIT: If I try to evaluate $I_>$ from joriki's answer, I note that $q$ is small and $x$ is large, so $qx$ can be anything, but $\epsilon=q\left(\sqrt{1+x^2-2x\cos\phi}-x\right) = q[-\cos\phi + \sin^2\phi/2x + O(x^{-2})]$ is small, so
\begin{align}
I_> &= q\int_1^\infty dx \int_0^{2\pi} \frac{d\phi}{2\pi}  U(qx +\epsilon) - U(qx) \\
&= \sum_{n=1}^\infty \frac{q^{n+1}}{n!}\int_1^\infty dx \, U^{(n)}(qx) \int_0^{2\pi} \frac{d\phi}{2\pi}   \left(\sqrt{1+x^2-2x\cos\phi}-x\right)^n \\
&= \frac{q^2}{2}\int_1^\infty dx \, U'(qx) \int_0^{2\pi} \frac{d\phi}{2\pi}   \left(\sqrt{1+x^2-2x\cos\phi}-x\right)\\ 
&+ \frac{q^3}{6}\int_1^\infty dx \, U''(qx) \int_0^{2\pi} \frac{d\phi}{2\pi}   \left(\sqrt{x^2+1-2x\cos\phi}-x\right)^2 + O(q^4)\\
&= \frac{q^2}{2}\int_1^\infty dx \,  U'(qx) f(x) + \frac{q^3}{6}\int_1^\infty dx \, U''(qx) [1-2x f(x)] + O(q^4)\\
f(x) &= \frac{(x+1)E[4x/(x+1)^2] + (x-1)E[-4x/(x-1)^2]}{\pi} - x = \frac{1}{4x} + O(x^{-3})\\
\end{align}
so, since $U'(q) = -\frac{2 q}{\sqrt\pi} + (1+2 q^2)e^{q^2} \erfc(q)$, the first $x$ integral diverges logarithmically and the next one is worse. Expanding in $\epsilon$ gives an expansion in $q$ which is divergent. Perhaps someone has an idea of how to proceed to get an expansion in $1/q$ for small $q$. 
 A: $\DeclareMathOperator{\erfc}{erfc}$Here's my answer so far; I'll try to complete it when I have the time, hopefully in a couple of days.
The problem is basically that you have difficult behaviour both for $q'\ll q$ and for $q'\gg q$. For $q'\ll q$, the integrand near $\phi=0$ behaves like $|q'-q|$, which isn't differentiable at $q'=q$. This problem could be avoided by expanding $U$ into a power series at $0$ instead of expanding in a power series in $q$, but this runs into the problem for $q'\gg q$: Here $U$ tends to a constant, but this isn't visible in a series expansion at $0$.
The solution, I believe, is to split the integral into two parts with $q'\lessgtr q$, much as one does in the Laplace expansion, in which the same argument occurs. So
\begin{align}
f(q)
&=
\int_0^\infty\mathrm dq'I(q,q')
\\
&=
\int_0^\infty\mathrm dq'\int_0^{2\pi} \frac{\mathrm d\phi}{2\pi}\left(U\left(\sqrt{q^2 + q'^2 -2qq'\cos\phi}\right)-U\left(q'\right)\right)
\\
&=
q\int_0^\infty\mathrm dx\int_0^{2\pi} \frac{\mathrm d\phi}{2\pi}\left(U\left(q\sqrt{1 + x^2 -2x\cos\phi}\right)-U\left(qx\right)\right)
\\
&=
\overbrace{q\int_0^1\mathrm dx\int_0^{2\pi} \frac{\mathrm d\phi}{2\pi}\left(U\left(q\sqrt{1 + x^2 -2x\cos\phi}\right)-U\left(qx\right)\right)}^{I_\lt}
\\
&+
\underbrace{q\int_1^\infty\mathrm dx\int_0^{2\pi} \frac{\mathrm d\phi}{2\pi}\left(U\left(q\sqrt{1 + x^2 -2x\cos\phi}\right)-U\left(qx\right)\right)}_{I_\gt}\;.
\end{align}
For $I_\lt$, we can expand $U(t)=t\mathrm e^{t^2}\erfc(t)=t-\frac2{\sqrt\pi}t^2+t^3+O\left(t^4\right)$ and integrate term by term:
\begin{align}
q^2\int_0^1\mathrm dx\int_0^{2\pi} \frac{\mathrm d\phi}{2\pi}\left(\sqrt{1 + x^2 -2x\cos\phi}\,-x\right)&\approx0.587336q^2
\end{align}
(Wolfram|Alpha quadrature),
\begin{align}
q^3\int_0^1\mathrm dx\int_0^{2\pi} \frac{\mathrm d\phi}{2\pi}\left(\left(1 + x^2 -2x\cos\phi\right)-x^2\right)&=q^3\;,
\end{align}
and so on, with the terms corresponding directly to powers of $q$.
I haven't solved the $I_\gt$ part yet, but I suspect that it has a non-analytic behaviour like $\exp\left(-\frac1{q^2}\right)$ at $q=0$, so it or a part of it needs to be expanded in $\frac1q$, so the overall result might be a sum of a series in $q$ and a series in $\frac1q$.
A: $\DeclareMathOperator{\erfc}{erfc}$A simple start
(which you may have already done):
Since
$\erfc(q)
\sim \dfrac{e^{-q^2}}{q\sqrt{\pi}}
$
for large $q$,
$U(q) 
=  q e^{q^2}\erfc q
\sim  q e^{q^2}\dfrac{e^{-q^2}}{q\sqrt{\pi}}
=  \dfrac{1}{\sqrt{\pi}}
$
so that
$U(q)$
is constant for large $q$.
For small $q$,
$\erfc(q)
= 1-\dfrac{2q}{\sqrt{\pi}}+O(q^3)
$
so that
for small $q$,
$\begin{array}\\
U(q) 
&=  q e^{q^2}\erfc q\\
&\sim  q e^{q^2}(1-\dfrac{2q}{\sqrt{\pi}}+O(q^3))\\
&\sim  q (1+q^2+O(q^4))(1-\dfrac{2q}{\sqrt{\pi}}+O(q^3))\\
&\sim  q (1-\dfrac{2q}{\sqrt{\pi}}+q^2+O(q^3))\\
\end{array}
$
By seeing the range of arguments
to $U$
in the integral,
and possibly using
more terms of the
asymptotic expansion
of $\erfc$,
you might get
some useful approximations.
