A novel method is proposed for flexibility analysis of nonconvex heat exchanger networks. In this method, the direction matrix is introduced to describe the deviation direction of the uncertain parameter. Then by searching the critical directions that restrict the process flexibility with the simulated annealing algorithm, the flexibility index can be obtained. Since the directions are no longer limited to the vertices, this method can well deal with the nonconvex problems. Moreover, for large-scale problems, a decoupling strategy is developed to enhance the efficiency. On the basis of this, the entire network is decomposed into several independent subnetworks, and the one with the lowest flexibility decides the global flexibility index. Two examples are studied, and the results well demonstrate the various aspects dealt with in this work.