During food freezing, controlling the crystallization is a key factor in quality. Most often, modelling of freezing is carried out assuming local thermodynamic equilibrium, considering the food product as an equivalent monophasic material and solving numerically the problem with cells much larger than the ice crystals. In the present work, a numerical model was developed at the scale of the ice crystals with cells of some micrometers. This way, the solute/water mutual diffusion and the conduction of latent heat release in the vicinity of each crystal is directly taken into account. The model was implemented with Matlab on a 2D mesh by considering in each cell the temperature (energy equation) and the volume fraction of water, ice and solutes (phase change and diffusion equations). The model allows to analyze and explain the influence of heterogeneous nucleation sites and surrounding temperature on the evolution of ice fraction and crystals size distribution.