Currently, RF is actively developing the Northern territories. Questions of studying the physical processes of icing are relevant, since climate conditions affect the surface of the objects under study (power lines, residential buildings, power plants, aircraft), human safety and ecology. In clouds, the appearance and movement of liquid droplets-particles is possible. When studying two-phase flows containing a suspension of aerosol particles (dispersed phase) in the carrier medium (dispersion medium) in the atmosphere, it is important to correctly evaluate the main parameters that define the system, and adequately describe the real process using a formulated mathematical model. This article is devoted to the development of a new iceFoam solver as part of the OpenFOAM v1912 package for modeling the icing process at a typical particle size of about 40 microns, which corresponds to Annex C of the AP-25 Aviation rules. The Euler-Lagrangian approach and finite volume method are used to describe the dynamics of liquid droplets. A modified liquid film model based on the shallow water theory is used as a thermodynamic model. The results of calculation for the case of flow around the cylinder and airfoil NACA 0012 using the URANS method and Spalart-Allmaras turbulence model are presented. In the calculation domain, two variants of grids are constructed: for an external gas-drop flow and for a liquid thin film with a thickness of one cell in height. Patterns of ice thickness distribution are given. When developing the source code using C++ language, the technology of inheritance was used, i.e. creating base and derived classes. As a result, a parallel iceFoam solver was developed for simulating the motion of liquid particles and the formation of ice on the bodies’ surface. For the calculation of one test case 8-32 computing cores were used on the ISP RAS HPC.