One of the most important requirements for the materials (zirconium alloys) used in the reactor active zone is low hydrogen absorptivity, since hydrogen-induced embrittlement may cause zirconium cladding damage. Depending on the hydrogen content and operation temperature, hydrogen may be present in zirconium alloys as a solid solution or as hydrides. Hydrides have the greatest embrittlement effect on alloys as they can initiate and enlarge cracks. The problem is to model the dynamics of the moving boundary of phase transition and to estimate the concentration distribution in the hydride and the solution. This paper presents a mathematical model of zirconium alloy hydrogenation taking into account the phase transition (hydride formation) and the iterative computational algorithm for solving the nonlinear boundary-value problem with free phase boundary based on implicit difference schemes. The study was carried out under state order to the Karelian Research Centre of the Russian Academy of Sciences (Institute of Applied Mathematical Research KarRC RAS).