The magnetohydrodynamic (MHD) activity during density limit disruptions in tokamaks is modelled numerically by three-dimensional resistive reduced MHD simulations with a simple transport model including radiation losses. The simulations reproduce experimentally observed phenomena such as the destabilization of MHD modes near the plasma edge during the early profile contraction phase, followed by growth of the m = 2/n = 1 mode to large amplitude, a sequence of minor disruptions and the major disruption. A new theoretical model is given for the major disruption, which takes place in two phases: (1) an internal relaxation flattens the temperature in the central part of the discharge and (2) the current profile broadens. The internal instability of the first phase has a mainly m = 1/n = 1 convection pattern, but, because of non-linear coupling to the large m = 2/n = 1 mode, the magnetic perturbation has a strong m = 3/n = 2 component. During the internal relaxation, the large amplitude 2/1, 1/1 and 3/2 perturbations break up the magnetic surfaces isolating the q almost-equal-to 1 region from the stochastic region around q = 2, and the magnetic field becomes stochastic in the entire q less-than-or-equal-to 2 region. In the second phase of the major disruption, MHD turbulence first develops on the stochasticized fields, resulting in current filamentation, initially in the central region where q less-than-or-equal-to 2. This leads to a broadening of the central current profile and a strong instability of the 2/1 mode. The disruption ends with rapid growth of the m > 2/n = 1 modes. The result is stochastic magnetic fields across the entire plasma and a large scale broadening of the current profile.