Coaxial high-speed gas jets are widely used to atomize paint liquids in industrial painting processes, such as using pneumatic or HVLP (High Volume Low Pressure) spray guns, where the liquid jet is broken up in a high-speed annular gas stream. Usually, a critical state with sonic speed does exist around the exit of the liquid jet. At the same time, a so-called shaping gas flow is used to form the spray pattern. Although some research work has been done to investigate the primary breakup of coaxial liquid-gas jets, prediction of full breakup processes is still an open research topic, since the simulation of a complete spray process using Volume of Fluid (VOF) method requires very fine local grids, which makes it unrealistic for practical applications. Conventionally, in order to perform spray painting simulation, namely the prediction of the droplet transport and the resulting film thickness distribution on a workpiece, empirical assumptions for the initial conditions of the droplet phase at the inlet plane are required, i.e., positions, velocity vectors and concentrations. Furthermore, droplet size distributions measured at locations further downstream, where droplet formation is fully developed, had to be applied. The purpose of the present study is to carry out numerical simulation of the atomization of viscous liquids, focusing on the combination of the liquid breakup and droplet transport processes. The numerical simulations were carried out using the commercial CFD code ANSYSFluent. The VOF-to-DPM model, namely coupled Volume of Fluid and Lagrangian particle tracking approaches, was used. Secondary breakup model was also applied. Effect of model parameters and grid size on the droplet distribution was studied. The intact liquid length close to the liquid nozzle and the formation of the liquid droplets in the regions of primary and fully developed breakup along the spray jet direction were analysed. Droplet size distributions downstream the nozzle were compared with experimental results.