Traditional methods, including direct solution methods based on Newton's method and indirect solution methods based on thermodynamic principles, are the mainstream methods used to solve the volume-temperature flash calculation (called NVT-flash), even though they suffer from drawbacks such as sensitivity to initial value and complexity of derivative calculations. A constrained backtracking search algorithm (CBSA), proposed in 2024, was the first and only metaheuristic algorithm to successfully tackle the NVT-flash problem, which overcomes shortcomings of traditional methods. Considering the advantages of metaheuristic algorithms, a constrained gray prediction evolutionary algorithm with a surrogate model based on quadratic interpolation (CGPE-QI) is proposed in this paper to deal with the NVT-flash problem. CGPE-QI considers total Helmholtz free energy as the objective function, moles vector, and volume of a single phase as variables. Constraints to solve the NVT-flash problem are addressed by using a direct search method and an exterior point method. Numerical experiments on two-phase equilibrium of pure substance and mixtures are carried out employing CGPE-QI. Experimental results are the same as those obtained by traditional methods, which confirms that CGPE-QI can effectively tackle the NVT-flash problem and possesses energy decay property. In particular, the results demonstrate that CGPE-QI is more competitive than CBSA in terms of convergence speed, stability, and calculation cost. CGPE-QI proposed in this paper is the second metaheuristic algorithm to successfully solve the NVT-flash problem, illustrating that metaheuristic algorithms have great potential in solving phase equilibrium calculation problems.