Recently, a number of studies have reported the formation of atmospheric superrotation in three-dimensional Venus GCMs. However, its mechanism has not been fully understood, which is partly due to the difficulty of model data analysis and possible numerical problems (Lebonnois et al 2012). Earlier studies with two-dimensional models by Matsuda (1980, 1982) and Yamamoto and Yoden (2013) provided useful insights, but to assess their applicability is not straightforward, since these studies are based on the Boussinesq approximation and use fixed uniform eddy diffusivity and Newtonian cooling rate. Here, I present the results of numerical simulations with a two-dimensional axi-symmetric model to overcome these difficulties. The model is based on the transformed Eulerian mean (TEM) equations supporting deep atmosphere over multiple scale heights, and it conserves total zonal angular momentum very well when surface friction is tuned off. The model does not use prescribed eddy diffusivity. Instead, effects of dynamical instabilities (convective, symmetric, barotropic-baroclinic instabilities) are parameterized. The model was initialized with the uniform solid-body rotation at the rate of the planetary rotation and run under radiative forcing similar to that in the recent Venus AFES simulations. A realistic superrotation emerged, and it reached a statistical steady state. It is revealed that the dynamical disturbances associated with dynamical instabilities can work like eddy diffusivity in the Gierasch mechanism to create a local maximum of the angular momentum within the atmosphere. The steady-state meridional circulation is multi-layered. The surface friction works to keep the angular momentum at low levels. Effects of thermal tides will be investigated elsewhere.