In this paper, a numerical approach is developed to investigate the evolution of fracture process zone (FPZ) during the complete fracture process in concrete structures by using stress intensity factor-superposition method. In this approach, the initial fracture toughness KICini, as an inherent material property, is introduced to form a crack propagation criterion for concrete. The developed numerical approach is then employed to analyze the complete fracture process of three series of notched concrete beams under three-point bending. It is found that the numerical results agree well with experimental ones published in literature through which the developed numerical approach, with an initial fracture toughness based crack propagation criterion, for fracture analysis of concrete is validated. The verified numerical approach is then utilized to simulate the complete fracture process of a series of concrete square plates with different sizes and/or initial crack length-to-depth ratio (a0/D). The effects of a0/D on evolution of FPZ length (aFPZ), especially after the FPZ is fully developed, are examined based on numerical analysis results. It is found that there are three different types of aFPZ variation with respect to a0/D, viz. (i) aFPZ keeps increasing after FPZ is fully developed. (ii) aFPZ turns to decrease from the peak value after FPZ is fully developed. (iii) FPZ is not fully developed. Finally, features of KR-curve for concrete are investigated based on the developed numerical method, and it is found that KR-curve of concrete is size-dependent when the effects of FPZ variation are taken into account.