In this work, we calculate the sub-leading power contributions to radiative leptonic D ->gamma lv decay. For the first time, we provide the analytic expressions of next-to-leading power contributions and the error estimation associated with the power expansion of O(Lambda(QCD)/m(c)). In our calculation, we adopt two different models of the D-meson distribution amplitudes phi(+)(D,I) and phi(+)(D,II). Within the framework of QCD factorization as well as the dispersion relation, we evaluate the soft contribution up to the next-to-leading logarithmic accuracy and also consider the higher-twist contribution from the two-particle and three-particle distribution amplitudes. Finally, we find that all the sub-leading power contributions are significant at lambda(D)(mu(0)) = 354 Mev, and the next-to-leading power contributions lead to 143% in phi(+)(D,I) and 120% in phi(+)(D,II)corrections to leading power vector form factors with E-gamma = 0.5 GeV. As the corrections from the higher-twist and local sub-leading power contributions are enhanced with increasing inverse moment, it is difficult to extract an appropriate inverse moment of the D-meson distribution amplitude. The predicted branching fractions are (1.88(-0.29)(+0.36)) x 10(-5) for phi(+)(D,I) and (2.31(-0.54)(+0.65)) x 10(-5) for phi(+)(D,II)