This study focuses on the simulation method for crystallization processes, that is the most important unit operations and is widely used in the chemical industry such as pigments, pharmaceuticals, and detergents, etc. The process includes the kinetics of nucleation, growth, agglomeration, and breakage. Because the models are very complex to become partial and ordinary differential, algebraic and integral equations(IPDAE), it is very difficult to get a numerical solution. And, it makes process design and analysis depend on experiments. In this study, the simulation is based on the method of characteristics to avoid the difficulties of numerical diffusion and stability that bother the previous methods in this area. In order to improve the accuracy and efficiency of the simulation, this study combines the method of characteristics with the adaptive mesh method, which adaptively adds and eliminates meshes as the solution varies. The proposed method can conserve crystal number and mass in the representation of agglomeration/breakage event and mesh adaptation. The proposed method is validated through the simulation of case studies.