In this talk we present an algorithm based on the Boundary Element Method for solving the Stokes flow around an arbitrarily shaped particle in flow. Assuming the particle is smaller than the Kolmogorov length scale of turbulence, the objective of the algorithm is to determine the drag force and torque experienced by the particle during its motion in a turbulent flow field. The algorithm is based on analytical expressions for the computation of singular integrals obtained in the discretisation of boundary integral equations for the Stokes problems with triangular or quadrilateral boundary elements. The algorithm is implemented for parallel execution while using linear interpolation of functions and constant boundary elements for fluxes. It is able to handle complex particle shapes with high accuracy and efficiency. First, we limit ourselves to superelliposid particle shapes and by carefully choosing flow conditions around the particle run tens of thousands of BEM simulations for different shapes and flow conditions. Using this simulation results we are able to create a neural network based surrogate model that is able to predict the drag and torque on the particle for a wide range of particle shapes and orientations with respect to the flow direction. Second, we study particle interactions for cases of higher volume fractions where hydrodynamic interactions between particles cannot be neglected. We show that by using the BEM solver we are able to accurately estimate drag and torque and simulate the motion of up to ten particles in proximity. Finally, we couple the developed BEM based surrogate model with the OpenFOAM CFD solver to simulate the motion of particles in a turbulent flow field. We study the motion of particles of different shapes and sizes in a human respiratory track model and show the influence of particle shape on deposition patterns.