This thesis presents the development and simulation of batch and continuous crystallization models. Especially, models are derived for simulating batch and continuous enantioselec- tive preferential crystallization processes in single and coupled crystallizers. Such processes are highly important in chemical and pharmaceutical industries. The effects of nucleation, growth, and fines dissolution phenomena on the crystal size distribution (CSD) are inves- tigated. For the first time continuous preferential crystallization is investigated and the effects of different seeding and operating strategies on the process are analyzed. To judge the quality of the process some goal functions are used, such as purity, productivity, yield and mean crystal size of the preferred enantiomer. The semi-discrete high resolution finite volume schemes (HR-FVS) and the discontinuous Galerkin (DG) finite element method are proposed for solving these models. The resulting systems of ordinary differential equa- tions (ODEs) are solved by using explicit and nonlinearly stable high order Runge- Kutta method. The schemes satisfy the total variation bounded (TVB) property which guarantees the positivity of the schemes, for example the non-negativity of CSD in the present case. The suggested methods have capabilities to capture sharp discontinuities and narrow peaks of the CSD. In DG-schemes, the accuracy of the method can be improved by introducing additional nodes in the same solution element and, thus, avoids the expansion of mesh stencils which is normally observed in high order finite volume schemes. For that reason, the method can be easily applied up to boundary cells without loosing accuracy. It was found that the proposed numerical schemes have the capability to solve the given models more efficiently and accurately. The results support process design and optimization.