the following problem $$\label{1.1}
\begin{cases}
-{\varepsilon}^2 \Delta u=\left(u-q-\frac{\kappa}{2\pi}\ln\frac{1}{{\varepsilon}}\right)_+^{p},  & \text{in}\;{\Omega},\\
u=0, &\text{on}\; \partial{\Omega},
\end{cases}$$ and gave exact asymptotic behavior and expansion of the least energy solution by estimating the upper bounds on the energy. The solutions for in [@SV] were obtained by finding a minimizer of the corresponding functional in a suitable function space, which can only give approximation to a single point non-vanishing vortex. This method is hard to obtain multiple non-vanishing solutions.

In this paper, we approximate stationary vortex solutions of Euler equations with multiple non-vanishing vorticity by stationary classical solutions.

Our main result concerning is the following:

\[nth1\]

Suppose that $\Omega\subset \mathbb{R}^2$ is a bounded simply-connected smooth domain. Let $v_n:
\partial\Omega\rightarrow \mathbb{R}$ be such that $v_n\in
L^s(\partial\Omega)$ for some $s>1$ satisfying $\int_{\partial\Omega}v_n=0$. Let $\kappa_i>0,~i=1,\cdots,m$. Then, for any given stable critical point $(x_1^*,\cdots,x_m^*)$ of Kirchhoff-Routh function $\mathcal{W}(x_1,\cdots,x_m)$ defined by , there exists ${\varepsilon}_0>0$, such that for each ${\varepsilon}\in
(0,{\varepsilon}_0)$, problem has a stationary solution $\mathbf{v}_{\varepsilon}$ with outward boundary flux given by $v_n$, such that its vorticities $\omega_{\varepsilon}$ satisfying $$supp (\omega_{\varepsilon})\subset\cup_{i=1}^m
B(x_{i,\,{\varepsilon}},C{\varepsilon})~~\text{for}~~x_{i,\,{\varepsilon}}\in\Omega,
~~i=1,\cdots,m$$ and as ${\varepsilon}\rightarrow0$ $$\int_\Omega \omega_{\varepsilon}\rightarrow \sum_{i=1}^m\kappa_i,$$

$$(x_{1,\,{\varepsilon}},\cdots,x_{m,\,{\varepsilon}})\rightarrow (x_1^*,\cdots,x_m^* ).$$

The simplest case, corresponding to a single point vortex ( $m=1$ ) was studied by Smets and 