We present a new software package, called QMix, for simulating the behavior of Superconductor / Insulator / Superconductor (SIS) mixers. The software uses a harmonic balance procedure to calculate the AC voltage across the SIS junction and multi-tone spectral domain analysis (MTSDA) to calculate the quasiparticle tunneling currents. This approach has two major advantages over other simulation techniques: (1) it can include an arbitrary number of higher-order harmonics, and(2) it can simulate the effect of multiple strong non-harmonic frequencies. This allows the QMix software to simulate a wide range of SIS mixer operation, including the effects of harmonics in the local-oscillator signal and gain saturation with respect to RF signal power.In this paper, we compare the simulated results from QMix to the measured performance of a 230 GHz SIS device, both to validate the software and to investigate the experimental data. To begin, we simulated the conversion gain of the mixer and we found excellent agreement with the experimental results. We were also able to recreate the broken photon step effect by adding higher-and lower-order harmonics to the LO signal. We then simulated a range of incident RF signal powers in order to calculate the gain saturation point. This is an important metric for SIS mixers because we require linear gain for reliable measurements and calibration. In the simulated results, we found both gain compression and gain expansion, which is consistent with other studies. Overall, these examples demonstrate that QMix is a powerful software package that will allow researchers to simulate the performance of SIS mixers, investigate experimental results and optimize mixer operation. We have made all of the QMix software open-source and we invite others to contribute to the project.