One of the most important requirements for the reactor's active zone materials (made of zirconium alloys) is low hydrogen absorptivity since hydrogen 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 form and enlarge cracks. The problem is to model the dynamics of the moving boundary of phase transition and to estimate the concentration distribution in hydride and in solution. This paper presents a mathematical model of zirconium alloy hydrogenation taking into account phase transition (hydride formation), the iterative computational algorithm for solving the nonlinear boundary-value problem with free phase boundary based on implicit difference schemes, and the results of computational experiments.