The multiple effect evaporator (MEE) system is a typical process for liquor concentration in the energy-intensive industries. Many approaches and commercial software have been widely used for MEE simulation and optimization. However, the existing work usually assumes simplified correlations for material thermodynamic properties and evaporator operations to avoid computational complexity, which may cause a large deviation in results. This Article addresses accurate and detailed unit operations and material properties to model the MEE systems and further heat integration for optimal energy recovery. To deal with the resulting complex mixed integer nonlinear programming (MINLP) problems of MEE systems, an efficient optimization strategy is developed, where the noncritical variables in the MINLP model are initialized as parameters and updated by solving a series of mixed integer linear programming (MILP) problems using a two-stage iterative procedure. An industrial scale problem for concentrating black liquor in a Chinese paper mill is carried out to demonstrate the validity and efficiency of the new approach. On the basis of the conditions of constant heat-transfer coefficients and stream boiling point rises assumed in the well-known commercial software WinGEMS, our method performs identically to WinGEMS in five distinct scenarios. Moreover, our method is more capable of solving industrial problems in practical situations including varying stream thermal properties and evaporator heat-transfer coefficients, achieving up to 25% of energy conservation in a real-world case.