Suspension plasma spraying process, a relatively new deposition technique in thermal spray coating, has been increasingly applied to deposit high-quality thermal barrier coatings using submicron particles. An accurate simulation of the process includes the development of a realistic model of the plasma both within and outside the torch. In this work, a three-dimensional time-dependent model has been developed to simulate the magnetohydrodynamic fields inside a DC plasma torch including arc fluctuations. The Reynolds stress model is used to simulate the time-dependent turbulent plasma flow. To investigate the effects of plasma arc fluctuation on the trajectory, temperature, and velocity of suspension droplets injected into the plasma jet, a two-way coupled Eulerian–Lagrangian method is employed. Submicron yttria-stabilized zirconia particles, suspended in ethanol, are modeled as multicomponent droplets. The Kelvin–Helmholtz Rayleigh–Taylor breakup model is used to simulate the droplet breakup. Particles are also tracked after the completion of suspension breakup and evaporation to obtain the in-flight particle conditions including the trajectory, size, velocity, and temperature. The arc attachment spots showed a good agreement with the experimental images. It was also shown that the properties of the particles are significantly affected by plasma arc fluctuations.