The well-known theory of absorption and fluorescence is briefly reviewed in a systematic manner for the Na D transitions. The resulting formalism is applied to simulation of Doppler-free saturation fluorescence spectra. With only one adjusting parameter, the nonradiative rate chosen to represent the time a thermal atom takes to move across the laser beams, the simulated Doppler-free spectra match the measured ones well for both D-1 and D-2 transitions over one-decade of excitation intensities. Relative to the weighted center of the six D-2 hyperfine transition lines, the frequencies of the dominant Doppler-free features have been determined from a simulated spectrum to within +/-0.1 MHz to be -651.4, 187.8, and 1068.0 MHz, respectively, for D-2a, crossover, and D-2b resonances. These features may be used as accurate frequency references for atmospheric spectroscopy. They are essential for the operation of the newly developed narrow-band Na fluorescence lidar for wind and temperature measurements in the mesopause region.