Controlling the grain size is an essential issue in the metal plastic working processes. Based on the consensus view that grain boundaries are impassable borders for dislocations and that the number of dislocations within a grain affects the stress builds up among the adjacent grains, altering grain size has a great influence on the post-forming mechanical properties of the deformed workpieces. In the present research, the deformation behaviors of 20CrMoH, a typical hardenability steel, are investigated and the microstructure evolution of grain growth experiencing heat load and grain refinement under deformation load is also uncovered. According to the experimental stress-strain curves and grain size variation data, a set of internal variable constitutive equations are established, in which the material constants are calibrated by a genetic optimization algorithm. Comparing the predicted curves with the experimental values, it can be concluded that the constitutive equations can accurately reproduce the thermal deformation behaviors and grain size variation of the steel. The correlation coefficient of the calculated and experimental values is 0.9861, and the average relative error approximates 4%. The unified constitutive equations are compiled into a commercial finite element software with the help of user-defined subroutine interfaces. The simulation results demonstrate that the maximum accuracy error of a hollow shaft size is about 10%, and the error in the average grain size is 12.1%. It is confirmed that the internal variable constitutive equations describing the evolution of grain size established in this paper are suitable for the numerical simulation of the bulk-forming process.