This article presents descriptions of theoretical models and numerical methods for simulating turbulent thermal plasma flow with nanopowder growth. Turbulence models must express turbulent and laminar states because both states co-exist with thermal plasmas showing large density variation and transport properties. Time-dependent 3D simulations are conducted based on Large Eddy Simulation using a dynamic Smagorinsky model. Results show significant difference depending on temporal and spatial discretization schemes and velocity-pressure coupling algorithms. Simulation results demonstrate that advanced numerical methods with high-order accuracy should be used for long and robust computations capturing steep gradients of nanopowder concentration and plasma temperature and 3D dynamic motions of multiscale vortices, which are turbulent features of thermal plasma flows with low Mach numbers. A thermal plasma jet generates a double-layer structure of inner high-temperature thick vortex rings and outer low-temperature thin vortex rings near the nozzle exit. Flowing downstream, these vortices interact, deform, and break up. Consequently, plasma transits to a complex thermal flow. The widely spreading distribution of multiscale vortices agrees with experimental observations, which are not simulated using conventional methods. Nanopowder is generated from material vapour by nucleation and condensation at interfacial regions between plasma and cold gas. Those regions include numerous vortices. Therefore, the vortices convey the nanopowder, producing a complex distribution of nanopowder. Simultaneously, the nanopowder diffuses and increases in size, decreasing in number by interparticle coagulation. Cross-correlation analysis suggests that a nanopowder distribution distant from a plasma jet can be controlled through temperature fluctuation control at the upstream plasma fringe.