We develop a new effective mathematical formulation and resulting universal computational algorithm capable of tackling various HF models in the framework of a unified approach. The scheme is not limited to any particular elasticity operator or crack propagation regime. Its basic assumptions are: (i) proper choice of independent and dependent variables (with the direct utilization of a new one – the reduced particle velocity), (ii) tracing the fracture front by use of the Stefan condition (speed equation), which can be integrated in closed form and provides an explicit relation between the crack propagation speed and the coefficients in the asymptotic expansion of the crack opening, (iii) proper regularization techniques, (iv) improved temporal approximation, (v) modular algorithm architecture. The application of the new dependent variable, the reduced particle velocity, instead of the usual fluid flow rate, facilitates the computation of the crack propagation speed from the local relation based on the speed equation. This way, we avoid numerical evaluation of the undetermined limit of the product of fracture aperture and pressure gradient at the crack tip (or alternatively the limit resulting from ratio of the fluid flow rate and the crack opening), which always poses a considerable computational challenge. As a result, the position of the crack front is accurately determined from an explicit formula derived from the speed equation. This approach leads to a robust numerical scheme. Its performance is demonstrated using classical examples of 1D models for hydraulic fracturing: PKN and KGD models under various fracture propagation regimes. Solution accuracy is verified against dedicated analytical benchmarks and other solutions available in the literature. The scheme can be directly extended to more general 2D and 3D cases.