Perturbation Theory and Thermodynamic Limit Suppose we have a classical Hamiltonian that can be divided into an “easy” part $H_0$ and a “difficult” part $\Delta H$ that depends on a parameter $g$:
\begin{equation}
H = H_0 + g \Delta H ~.
\end{equation}
The partition function can be written as 
\begin{equation}
Z=\sum_\text{states}e^{-\beta H_0}e^{-g\beta \Delta H} ~,
\end{equation}
because it's a classical system so $H_0$ and $\Delta H$ commute. The partition function of the easy part,
\begin{equation}
Z_0=\sum_\text{states} e^{-\beta H_0} ~,
\end{equation}
is supposed to be easy to compute.
A common way to solve this problem is to work perturbatively on $g$. We find
\begin{equation}
\begin{aligned}
Z &= \sum_\text{states} \Big( e^{-\beta H_0}\Big\{ 1-g\beta \Delta H +\frac{1}{2}\big(g\beta \Delta H\big)^2 +\cdots\Big\} \Big) \\
&= \sum_\text{states} e^{-\beta H_0} \cdot \Big( 1 - g\beta\langle \Delta H \rangle_0 + \frac{g^2\beta^2}{2}\langle\Delta H^2\rangle_0 + \cdots \Big)
\end{aligned}
\end{equation}
Where $\langle A \rangle_0 \equiv \frac{1}{Z_0} \sum_\text{states} A \, e^{-\beta H_0}$ is the mean value of $A$ using $H_0$.
The problem is that $H_0 \propto N$ because it is an extensive variable, so $e^{-\beta H_0}\, \propto \ e^{N}$ while $\langle \Delta H \rangle_0$ only scales as $N$. When calculating the Helmholtz energy in the thermodynamic limit,
\begin{equation}
\begin{aligned}
f&= \lim_{N\rightarrow\infty}-\frac{kT}{N}\ln(Z) \\
&= \lim_{N\rightarrow\infty}-\frac{kT}{N}\ln(Z_0)-\frac{kT}{N}\ln\Big( 1 - g\beta\langle \Delta H \rangle_0 + \frac{g^2\beta^2}{2}\langle\Delta H^2\rangle_0 \Big) ~.
\end{aligned}
\end{equation}
The first term scales as $\frac{1}{N}\ln(e^{N})$ so it is finite, but the second term scales as $\frac{1}{N}\ln(1 + N + N^2)$, so it should vanish.
How do we use perturbation theory in the thermodynamic limit then? 
 A: As far as I know, in many-body condensed-matter theory, a simplistic finite-order Taylor expansion is (almost) never used, unless for pedagogical reasons. In fact, one always needs to “organize the series expansion more cleverly” [as higgsss commented] in form of infinite perturbative resummations or asymptotic expansions.
A simple Taylor-expansion, as in $ e^{-x} = 1 - x + \frac{1}{2}x^2 + \mathcal{O}(x^3) $, works only for some single-particle problems in physics. 
In words of Ref. [1] (§5.1, p. 195):

The moral to be taken ... is that perturbative expansions should not
  be confused with rigorous Taylor expansions. Rather they represent
  asymptotic expansions, in the sense that, for weaker and weaker coupling, a partial resummation of the perturbation series leads to
  an ever more precise approximation to the exact result. For weak
  enough coupling the distinction between Taylor expansion and
  asymptotic expansion becomes academic (at least for physicists).
  However, for intermediate or strong coupling theories, the asymptotic
  character of perturbation theory must be kept in mind.

Consult §5.1 of Ref. [1] for an extensive discussion.
In some contexts, eg. in the perturbative renormalization group, one actually “re-exponentiates” the finite-order perturbative expansion; for instance,
$$
Z = \sum_{\text{states}} e^{-\beta (E_0 + g \Delta E)}
= Z_0 \, \langle e^{-\beta \, g\Delta E} \rangle_0
\approx Z_0 \, \langle 1 - \beta g \Delta E \rangle_0
\approx Z_0 \, e^{- \beta g \langle \Delta E \rangle} ~,
$$
so that the problem observed in the question does not appear at all; see, for instance, §8.2 of Ref. [1].

[1] Altland, A., and B. D. Simons. “Condensed matter field theory”. CUP (2010) [wcat].
A: As you noticed, the first thing to do for the thermodynamic limit is to consider quantities which have a chance of making sense in this limit. The partition function $Z$ is not good for that, but the limit of $\frac{1}{N}\log Z$ (related to the free energy density or the pressure) is indeed the right quantity. Correlation functions also work. Using your notations, one has
$$
Z=Z_0 \langle e^{-g\beta\Delta H}\rangle_0
$$
so
$$
\log Z=\log Z_0+\log \langle e^{-g\beta\Delta H}\rangle_0
=\log Z_0+\sum_{m=1}^{\infty}\frac{(-g\beta)^m}{m!}
\langle \Delta H,\ldots,\Delta H\rangle_0^{\rm T}
$$
where $\langle \Delta H,\ldots,\Delta H\rangle_0^{\rm T}$ is the $m$-th joint cumulant of the random variable $\Delta H$ with itself (or simply the $m$-th cumulant of $\Delta H$). The underlying probability measure corresponds to expectations $\langle\cdots\rangle_0$. If you expanded the logarithm as one of the comments suggested, you would have seen, e.g., that the quadratic term would involves the variance
$$
\langle \Delta H,\Delta H\rangle_0^{\rm T}=\langle (\Delta H)^2\rangle_0-
\langle \Delta H\rangle_0^2\ .
$$
One may naively think that, since $\Delta H$ is intensive, the cumulant
$$
\langle \Delta H,\ldots,\Delta H\rangle_0^{\rm T}
$$
scales like $N^m$, but this is not true. It scales like $N$ for any $m$. 
There are several ingredients coming together for this property. The first is that a joint cumulant $\langle X_1,\ldots,X_m\rangle^{\rm T}$ is always zero if the random variables $X_1,\ldots,X_m$ can be arranged in two groups which are statistically independent. This forces a "connectedness" property. The second is that a decent $\Delta H$ should have a built-in decay property in the size and spatial extent of the support of the interaction. More concretely, take the Ising model with possibly long-range interactions. Then $H_0=0$ (if you want a less trivial $H_0$ turn on a magnetic field), and
$$
\Delta H=\sum_{x,y}J_{x,y} \sigma_x\sigma_y
$$
for some translation invariant couplings that satisfy a decay condition when $x,y$ get far apart.
In particular, one wants a condition like
$$
\sum_y |J_{x,y}|<\infty
$$
for $x$ fixed, otherwise $\Delta H$ would not be "extensive".
Inserting the definition of $\Delta H$ in the cumulant, you see that you have to sum over a configuration of $m$ edges thrown into the lattice. They must form a connected graph. One basically has to perform a sum
$$
\sum_{x_1,y_1}\cdots\sum_{x_m,y_m} 
J_{x_1,y_1}\cdots J_{x_m,y_m}
\langle \sigma_{x_1}\sigma_{y_1},\ldots,\sigma_{x_m}\sigma_{y_m}
\rangle_0^{\rm T}\ .
$$
This results in an overall factor of $N$ to pick the location of $x_1$ say.Then by translation invariance you can assume that this point is pinned at the origin say. But you have to sum over the $2m-1$ remaining points. This is where the connectedness and summability condition on the $J$'s will save you. 
This is a simple case of a very general method called a cluster expansion.
To learn more about this, have a look a the notes I wrote


*

*"Notes on the cluster expansion for the polymer gas, a.k.a., the Mayer expansion."

*"Notes on the Brydges-Kennedy-Abdesselam-Rivasseau forest interpolation formula" needed for Lemma 1 in the previous note and for more sophisticated versions of the cluster expansion.


This was from a graduate course I taught a while ago.
A: I think the problem arises because you assume that $\langle \Delta H^2\rangle_0$ is in the same order as $\langle \Delta H\rangle_0^2$, i.e in $O(N^2)$. In fact, for any canonical ensemble the energy fluctuation,
$$\frac {\langle H^2\rangle}{\langle H\rangle^2}\propto \frac{1}{N},$$
(see, for example in here). Thus, following your perturbation case, $\langle \Delta H^2\rangle_0$ should be in $O(N)$; you can check this at the higher orders (notice that this has some connections to Abdelmalek's answer).
Therefore for any $\left|\frac{g\Delta H}{H_0}\right|\ll1$ the partition function should converge .
A: Every answer was really useful so I will try to put it all together in one answer
As Abdelmalek Abdesselam said, the partition function can we written as
\begin{equation} \tag{1}\label{part}
    Z=\sum_{\{\sigma\}}e^{-\beta H_0}e^{-g\beta H'}=Z_0 \langle e^{-g\beta H'}\rangle_0
\end{equation}
Then, the free energy is
\begin{equation}
-\beta F = \log(Z)=\log(Z_0)+\log\langle e^{-g\beta H'}\rangle_0
\end{equation}
The last term can be expanded in Taylor series as
\begin{equation}
    \log\langle e^{-g\beta H'}\rangle_0=\log\Big(1-g\beta \langle H'\rangle_0 + \frac{g^2 \beta ^2}{2}\langle H'^2\rangle_0+...\Big)
\end{equation}
If now we expand the $\log$ in powers of $g$ we get
\begin{equation}
     \log\langle e^{-g\beta H'}\rangle_0=-g\beta\langle H'\rangle_0+\frac{g^2 \beta^2 }{2}\Big[ \langle H'^2 \rangle_0-\langle H' \rangle_0 ^2 \Big]+...=\sum_{n=1}^{\infty}\frac{(-g\beta)^n}{n!}C_{H'}^{\ (n)}
\end{equation}
Where $C_{H'}^{\ (n)}$ is the n-th momentum of the aleatory variable $H'$. The first one is the mean value, the second one is the variance $\langle \Delta H'^2 \rangle_0$ and so on. The free energy is then
\begin{equation} \tag{2}\label{ener}
   F = F_0+g\langle H'\rangle_0 - \frac{\beta g^2}{2}\langle \Delta H'^2 \rangle_0 + ...
\end{equation}
In equation (\ref{part}) one can easily see that $Z\ \alpha\ e^N$. So $F$ is clearly an extensive quantity. The question is if this holds at every order in perturbation theory. The answer is yes, because every momentum $C_{H'}^{\ (n)} $ is extensive so equation (\ref{ener}) holds to every order. The proof is fairly simple. If one notices
\begin{equation}
  \frac{1}{\beta^(n-1)} \frac{\partial^n F}{\partial g^n}|_{g=0}=C_{H'}^{\ (n)}
\end{equation}
because we've already proved that $F$ is extensive, derivatives respect to $g$ (which is intensive) will remain extensive, so every moment $C_{H'}^{\ (n)}$ is extensive too. 
This result make look familiar to those who have taken a course in statistical mechanics. As donnydm pointed out, recall that in the canonical ensemble the mean energy is obtained as
\begin{equation}
U=-\frac{\partial}{\partial \beta} \log(Z)
\end{equation}
and that the specific heat (an extensive quantity) is 
\begin{equation}
C_V=\frac{\partial U}{\partial T}=-\beta^2 \frac{\partial}{\partial \beta} U= \beta^2 \frac{\partial^2}{\partial \beta^2} \log(Z)
\end{equation}
The derivation can be performed explicitly and the answer is
\begin{equation}
C_V=\beta^2 \Big[\langle H^2\rangle -\langle H \rangle ^2 \Big]=\beta^2 \langle \Delta H^2\rangle
\end{equation}
The specific heat is proportional to the variance of the energy. This is just an example of the more general result that may be familiar to everybody.
