  The aim of this paper is to study the multiplicity of solutions for the
following Kirchhoff type elliptic systems
  \begin{eqnarray*}
  \left\{ \arraycolsep=1.5pt
  \begin{array}{ll}
  -m\left(\sum^k_{j=1}\|u_j\|^2\right)\Delta
u_i=\frac{f_i(x,u_1,\ldots,u_k)}{|x|^\beta}+\varepsilon h_i(x),\ \ & \mbox{in}\
\ \Omega, \ \ i=1,\ldots,k ,\\[2mm]
  u_1=u_2=\cdots=u_k=0,\ \ & \mbox{on}\ \ \partial\Omega,
  \end{array}
  \right.
  \end{eqnarray*} where $\Omega$ is a bounded domain in $\mathbb{R}^2$
containing the origin with smooth boundary, $\beta\in [0,2)$, $m$ is a
Kirchhoff type function, $\|u_j\|^2=\int_\Omega|\nabla u_j|^2dx$, $f_i$ behaves
like $e^{\beta s^2}$ when $|s|\rightarrow \infty$ for some $\beta>0$, and there
is $C^1$ function $F: \Omega\times\mathbb{R}^k\to \mathbb{R}$ such that
$\left(\frac{\partial F}{\partial u_1},\ldots,\frac{\partial F}{\partial
u_k}\right)=\left(f_1,\ldots,f_k\right)$, $h_i\in
\left(\big(H^1_0(\Omega)\big)^*,\|\cdot\|_*\right)$. We establish sufficient
conditions for the multiplicity of solutions of the above system by using
variational methods with a suitable singular Trudinger-Moser inequality when
$\varepsilon>0$ is small.
