Crystallization phenomena in stirred reactors are influenced by local hydrodynamic conditions and these must be taken into account for successful process scale-up and optimization. This article is the first of two presenting the application of the Smoothed Particle Hydrodynamics (SPH) method to modelling of batch crystallisation in stirred tanks. The benefits of the Lagrangian meshfree methods were discussed and the SPH method was proposed as an efficient method for the rapid prediction of the global mean flow in stirred reactors. Various aspects of the simulation results were discussed such as quality of the fluid prediction, computational requirements and availability of the crystal size distribution without reconstruction. It has been shown that the computational requirements and accuracy of the fluid flow prediction can be controlled through the particle size and that the particles with a radius of 1.50-1.75 mm in the reactor of 2.65 I provide a good balance between the quality of the prediction and the computational requirements. The developed Smoothed Particle Hydrodynamics model was applied to a numerical solution of coupled computational fluid dynamics and discretised population balance equations to model a batch crystallization process. Due to the specific formulation of the SPH equations the resulting ODE system is solved using the weighted contributions rather than numerically by solving a linear system of equations. Therefore, a large number of additional transport equations resulting from the discretisation of the population balance leads to only a minor increase in computational requirements (around 60% for 200 equations). The non-idealities in the reactor result in non-uniform mixing and contribute to the dispersion of the CSD. The effect of the hydrodynamics on the local temperature/supersaturation and the resulting crystal size distribution was captured and compared to the ideal mixing case. The developed SPH method served as a basis for the Part II of the series where a methodology for solution of population balance equations using high-resolution finite volume schemes or method of characteristics computed in parallel and independently from the Navier-Stokes equations is presented. (C) 2015 Elsevier Ltd. All rights reserved.