A numerical framework based on the finite-difference time-domain method is proposed for the rigorous study of electro-optically tunable terahertz devices based on liquid crystals. The formulation accounts for both the liquid-crystal full-tensor anisotropy and the dispersion of its complex refractive indices, which is described via modified Lorentzian terms. Experimentally characterized liquid-crystalline materials in the terahertz spectrum are fitted and modeled in benchmark examples, directly compared with reference analytical or semi-analytical solutions. In addition, the efficiency of broadband time-domain modeling of the proposed technique is also demonstrated by accurately reproducing time-domain spectroscopy measurements.