Mathematical models have been studied and solved for evaluating an optimized process configuration for the energy intensive black liquor concentrating Kraft recovery process in paper mills. In the present study, a heptads effect evaporator system is considered and modeled first on the basis of three possible flow directions of black liquor feed and heating steam, i.e. backward, forward or mixed feed. Further, live steam split, liquor feed split, feed preheating and a hybrid of these energy saving schemes are coupled with the basic backward, forward and mixed feed arrangements. The systematically evaluated material and heat balance equations evolve into main model equations that are then represented in matrix forms, thereby, to generalize the models mathematically. The advantage of the studied models are their simplicity in linear representation of equations in matrix form and ease of numerical solution. The proposed mathematical models are iteratively solved using different numerical techniques: Gauss-Jordan, Gauss-elimination, Gauss–Seidel, Jacobi, successive over-relaxation and interior-point methods. The simulation results indicate that amongst the 15 simulated models, backward feed with the added arrangements of feed split, steam split and feed preheating showed the best steam economy. The studied models can be applied and easily extended to solve problems with different operating conditions as well once the liquor, steam and other evaporator effects parameters are known.