Direct numerical simulations in fully three-dimensional geometry have been performed in order to investigate effects of both the parallel flow to the magnetic field lines discarded in the reduced equations and the parallel thermal conductivity on the nonlinear evolution of interchange instability in a low- magnetohydrodynamic equilibrium in the inward-shifted vacuum configuration of the Large Helical Device. As the parallel thermal conductivity becomes larger, the perpendicular velocity of the linear eigenfunction becomes relatively smaller than the parallel velocity, and in the nonlinear phase, the energy of the perpendicular velocity is selectively damped by the viscosity. Consequently, the plasma is finally saturated with flattened pressure profile, finite parallel velocity, and fairly good magnetic flux surfaces.