Difference frequency generation (DFG) is an essential phenomenon for generating far-infrared (FIR) and broadband terahertz (THz) waves, which cannot be generated by either alternating current flow-driven antennas for microwave source or pumped materials for higher frequency lights. Metallic nonlinear metasurfaces (NLMSs) as a new platform for generating THz sources were intensively developed. In this work, a nondepleted time-domain model derived from the Maxwell-hydrodynamic model is proposed where linear and nonlinear responses in metallic structures under an intense excitation are separated. By using the proposed model, only linear and low-frequency DFG signals will be considered, which makes this model suitable for unconditionally stable time-domain algorithms. A finite-difference time-domain (FDTD) method is employed to solve a set of multiphysics equations based on the proposed model. The stability of the algorithm is apparent. We validate our method by simulating THz generations from metallic NLMSs and benchmark the approach with previous numerical studies as well as our recent experimental results. In addition, the implementation of solving the DFG response can be remarkably accelerated by using the unconditionally stable algorithm.