Current simulations of swirling flows in pipes are limited to relatively low Reynolds number flows (Reβ<β6000); however, the characteristic Reynolds number is much higher (Reβ>β20,000) in most of engineering applications. To address this difficulty, this paper presents a numerical simulation algorithm of the dynamics of incompressible, inviscid-limit, axisymmetric swirling flows in a pipe, including the vortex breakdown process. It is based on an explicit, first-order difference scheme in time and an upwind, second-order difference scheme in space for the time integration of the circulation and azimuthal vorticity. A second-order Poisson equation solver for the spatial integration of the stream function in terms of azimuthal vorticity is used. In addition, when reversed flow zones appear, an averaging step of properties is applied at designated time steps. This adds slight artificial viscosity to the algorithm and prevents growth of localized high-frequency numerical noise inside the breakdown zone that is related to the expected singularity that must appear in any flow simulation based on the Euler equations. Mesh refinement studies show agreement of computations for various mesh sizes. Computed examples of flow dynamics demonstrate agreement with linear and nonlinear stability theories of vortex flows in a finite-length pipe. Agreement is also found with theoretically predicted steady axisymmetric breakdown states in a pipe as flow evolves to a time-asymptotic state. These findings indicate that the present algorithm provides an accurate prediction of the inviscid-limit, axisymmetric breakdown process. Also, the numerical results support the theoretical predictions and shed light on vortex dynamics at high Re.