Compositional simulation is challenging, because of highly nonlinear couplings between multi-component flow in porous media with thermodynamic phase behavior. The coupled nonlinear system is commonly solved by the fully-implicit scheme. Various compositional formulations have been proposed. However, severe convergence issues of Newton solvers can arise under the conventional formulations. Crossing phase boundaries produces kinks in discretized equations, and subsequently causing oscillations or even divergence of Newton iterations.The objective of this work is to develop a smooth formulation that removes all the property switches and discontinuities associated with phase changes. We show that it can be very difficult and costly to smoothing the conservation equations directly. Therefore, we first reformulate the coupled system, so that the discontinuities are transferred to the phase equilibrium model. In this way a single and concise non-smooth equation is achieved and then a smoothing approximation can be made. The new formulation with a smoothing parameter provides smooth transitions of variables across all the phase regimes. In addition, we employ a continuation method where the solution progressively evolves toward the target system.We evaluate the efficiency of the new smooth formulation and the continuation method using several complex problems. Comparing to the standard natural formulation, the developed formulation and method exhibit superior nonlinear convergence behaviors. The continuation method leads to smooth and stable iterative performance, with a negligible impact on solution accuracy. Moreover, it works robustly for a wide range of flow conditions without parameter tuning.