The magnetohydrodynamic (MHD) activity during density limit disruptions is modeled numerically by three-dimensional resistive reduced MHD simulations with a simple transport model and 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: first, an internal relaxation flattens the temperature in the central part and subsequently, the current profile broadens. The internal instability of the first phase has a mainly m=1/n=l displacement, but because of nonlinear coupling to the large m=2/n=1 mode, the magnetic perturbation has a strong m=3/n=2 component. At a late stage of the internal relaxation, the large amplitude 2/1, 1/1, and 3/2 perturbations give rise to stochastic magnetic fields in the q=l island, and at its end, the magnetic field is stochastic in the region q 2. During the second phase, MHD turbulence develops on the stochasticized fields, resulting in filamentation and broadening of the current density, first in the central region. The disruption ends with a rapid instability of the m2/n=1 modes. The result is stochastic magnetic fields across the entire plasma and a large-scale broadening of the current profile.