We perform two-dimensional numerical hydrodynamical calculations of thermal instability in cooling flows in a Cartesian coordinate system, in order to investigate the evolution of the infailing perturbation in the non-linear stage. The formation of vortices in the infailing perturbation is found. This is a new, non-linear phenomenon with important implications for cooling-flow structure and its evolution. Our numerical results are divided into two cases as follows. The first case is that the vortices are formed in both side parts of the perturbation by the non-linear effect. In this case, the amplitude of the density perturbation decreases and, as a result, it is prevented from radiative cooling. The second case is that the perturbation cools before the vortex is formed and, as a result, the low-temperature and high-density cold gas sheet is formed. The cold gas sheet is elongated in the direction of gravitational force. In this case, the vortices are also formed in the residual part of the perturbation after the cold gas sheet is formed. We give the criterion that the perturbation cools before the vortex is formed. Very large values of initial density contrast of perturbations are required for cooling to win over the vortex formation in real cooling flows. Then the apparently pervasive and extensive mass deposition occurring in cooling flows remains unexplained. We suggest that, if there are many perturbations which can evolve to the non-linear stage in a cooling flow, the cooling flow must be very turbulent because of the vortex formation and its development. If there is a weak magnetic field which freezes in the cooling flow, the magnetic field is very tangled by the vortex motion and then the effective mean free path and the conductivity would be much smaller than the Spitzer value.