When considering the time-variant characteristics of atmospheric turbulence in the numerical simulation research, it is very essential to generate the phase screens that meet the requirement. Meanwhile, in the numerical simulation of light propagation along an inhomogeneous turbulent path, the reasonable position of phase screens is also very important. In this article, an iterative algorithm of phase screen position based on the multilayer model of phase screens is given. And the phase screens are also continuously interpolated in order to meet the need of dynamic movement. Moreover, an example of light propagation combining the time-variant phase screens in two-layer system is provided under condition of the Hufnagel-Valley5/7. Thereafter, the characteristics of the optical wavefront and intensity changed with time are analyzed. The results show that the statistical properties of the generated phase screen coincides well with theoretical model of atmospheric turbulence in the low and high frequency region which is mainly concerned. The variance of Zernike coefficients of the optical wavefront in numerical simulation agrees with theoretical values basically. Although there is certain error existence in the probability density distribution of intensity fluctuation between simulation results and theoretical analyses, the trend of it is consistent with the lognormal distribution. Thus, the feasibility and validity of the above mentioned iterative algorithm is vindicated.