A significant percentage of computational time in compositional simulations is spent performing flash calculations to determine the equilibrium compositions of hydrocarbon phases in situ. Flash calculations must be done at each time step for each grid block; thus billions of such calculations are possible. It would be very important to reduce the computational time of flash calculations significantly so that more grid blocks or components may be used. In this dissertation, three different methods are developed that yield fast, robust and accurate phase behavior calculations useful for compositional simulation and other applications. The first approach is to express the mixing rule in equations-of-state (EOS) so that a flash calculation is at most a function of six variables, often referred to as reduced parameters, regardless of the number of pseudocomponents. This is done without sacrificing accuracy and with improved robustness compared with the conventional method. This approach is extended for flash calculations with three or more phases. The reduced method is also derived for use in stability analysis, yielding significant speedup. The second approach improves flash calculations when K-values are assumed constant. We developed a new continuous objective function with improved linearity and specified a small window in which the equilibrium compositions must lie. The calculation speed and robustness of the constant K-value flash are significantly improved. This new approach replaces the Rachford-Rice procedure that is embedded in the conventional flash calculations. In the last approach, a limited compositional model for ternary systems is developed using a novel transformation method. In this method, all tie lines in ternary systems are first transformed to a new compositional space where all tie lines are made parallel. The binodal curves in the transformed space are regressed with any accurate function. Equilibrium phase behavior calculations are then done in this transformed space non-iteratively. The compositions in the transformed space are translated back to the actual compositional space. The new method is very fast and robust because no iteration is required and thus always converges even at the critical point because it is a direct method. The implementation of some of these approaches into compositional simulators, for example UTCOMP or GPAS, shows that they are faster than conventional flash calculations, without sacrificing simulation accuracy. For example, the implementation of the transformation method into UTCOMP shows that the new method is more than ten times faster than conventional flash calculations.