With conventional reservoirs being depleted, unconventional reservoirs play a significant role in worldwide energy supply. Even though the exploitation of unconventional resources has achieved unprecedentedly great success, flow mechanisms and the underneath phase behavior are not clearly illustrated. Recent studies indicate that the confined space in the nanoscale and the considerable capillary pressure result in a significant deviation of fluid properties from their bulk properties. Failure to incorporate these effects will lead to inaccurate estimation of reservoir performance. As a consequence, the confinement effect and capillary pressure should be taken into account to accurately characterize the phase behavior in unconventional reservoirs. In this study, the capillary effect is incorporated into phase equilibrium calculations at constant moles, volume, and temperature, which is believed to have advantages over phase equilibrium calculations at constant moles, pressure, and temperature in many cases. We develop evolution equations for moles and volume using the laws of thermodynamics and Onsager’s principle. A thermodynamically stable numerical algorithm is constructed, which is the extension of our previous work, by introducing a generalized Onsager coefficient matrix. The new algorithm exhibits better performance in efficiency and accuracy by simultaneously solving the evolutionary equations. At the equilibrium state, the pressure inequality resulting from the capillary pressure is spontaneously satisfied. In addition, by properly selecting the initial approximation for phase equilibrium calculations, we obtain a physically meaningful transition process with uniform phase identification. A number of numerical examples is presented to demonstrate the robust performance of the proposed model. It is found that capillary pressure affects phase composition and distribution. The bubble point pressure is suppressed under capillary effect, while the dew point pressure exhibits complex behavior: increasing in the upper branch and decreasing in the lower branch of the dew point curve. In the limit of zero capillarity, our simulation results are consistent with the conventional isothermal–isochoric flash calculation.