Multi-stage alternating injection acid fracturing is commonly employed in the stimulation of tight carbonate reservoirs to enhance differential etching along the fracture surfaces and improve the conductivity of acid-etched fractures. The numerical simulation technique serves as an effective tool for optimizing the operational parameters of such treatments, significantly contributing to the enhancement of post-fracturing productivity and long-term production stability. However, existing numerical simulation approaches for multi-stage alternating injection acid fracturing often neglect acid-rock reactions or adopt simplified equivalent viscosity methods, which result in considerable deviations between simulation results and actual field observations. To address this issue, this study developed a mathematical model for multi-stage alternating injection acid fracturing based on the Volume of Fluid (VOF) method. This model incorporated both the interface tracking between reactive and non-reactive fluids, and also the acid–rock reaction. The governing equations of the mathematical model were discretized using the finite difference method, and the resulting numerical model was solved through computer programming. The accuracy of the model in capturing viscous fingering behavior and acid-etching profiles was verified by comparing the simulation results with the experimental data and analytical solutions. Based on this validated model, simulations were conducted to investigate the flow and reaction behavior of acid under different numbers of alternating injection stages, as well as the evolution of viscous fingering patterns and changes in etched fracture width. To comprehensively evaluate the effectiveness of differential etching, a viscous fingering index was introduced, which was accounted for acid penetration distance, the number of fingering branches, and the area covered by the viscous fingering. Simulation results demonstrate that under typical fracture widths and alternating injection conditions, low-viscosity acid gradually forms preferential flow channels in the fracture due to the viscosity contrast, which is the manifestation of the fingering phenomenon. As the number of alternating stages increases, the competitive development and mergence between the adjacent fingers happens. The effective acid penetration distance continues to increase with the number of alternating injection stages. However, when it is beyond a certain critical stage number, the growth rate of acid penetration distance slows, and further increasing of the alternating injection stages primarily only enhances the acid-etched width within the existing viscous fingering regions. Therefore, for a given fracture geometry and acid system, there is an optimal range of alternating stages, which simultaneously maximizes differential etching and the acid penetration distance. This model provides an effective simulation tool for the optimization of multi-stage alternating injection acid fracturing and offers theoretical guidance for the design of field treatment.