The present work is concerned with the development of mathematical models for the computer simulation of the dynamic operation of various industrial free-radical polymerization processes. For the development of the mathematical models, a generalized comprehensive kinetic mechanism was utilized for the description of the polymerization kinetics and the derivation of the population balance equations (PBEs) for the ‘live’ and ‘dead’ polymer chains. The correct solution of the PBEs can lead to the accurate calculation of the molecular characteristics of the polymers (e.g., molecular weight distribution (MWD), joint molecular weight – long chain branching distribution (MW-LCBD), etc.) that are directly related to the most significant mechanical and rheological characteristics of the polymer products. Much attention was given in the description of the diffusion phenomena presenting appropriate mathematical terms that quantitatively describe how these phenomena affect the polymer production and the molecular properties of the polymer produced. Furthermore, mathematical models were developed in order to describe the phase equilibrium calculations, the complex polymerization kinetic scheme that is produced by the use of bifunctional initiators and the dynamic behavior of industrial batch suspension polymerization reactors. The predictive capabilities of the proposed models are demonstrated by a direct comparison of model predictions and experimental data provided by the open literature and the polymer industry. A model-based study of the molecular weight distribution (MWD) of linear polymers and the combined molecular weight – long chain branching distribution (MW-LCBD) of branched polymers, using a population balance approach, is presented. In order to reduce the infinite system of differential PBEs, a lower system of differential equations can be obtained using the method of moments. This method is based on the statistical representation of the average molecular properties of the polymer (i.e., number and weight average molecular weights) in terms of the leading moments of the total number of chain length distributions (TNCLDs), for the ‘live’ and ‘dead’ polymer chains. To circumvent the well-known problem of closure of the ‘higher order’ moments, appearing in free-radical branched polymerization systems, the method of ‘bulk’ moments was employed. However, the method of ‘bulk’ moments cannot be used in processes where the reaction of random scission is present and cannot provide information on the full molecular weight distribution of linear or branched polymers. The orthogonal collocation on finite elements method (OCFE) and the fixed pivot technique (FPT) were used successfully in the calculation of the MWD of linear polymers. At first, the predictive capabilities of the methods were tested against experimental data on MWD. Then, the methods were compared with each other in terms of the accuracy and the computational effort spent for the calculations. Simulation results show that for linear polymers, the OCFE and FP methods result in similar molecular weight distributions for both unimodal and bimodal MWDs. However, when the quasi steady state approximation (QSSA) is used for all the radical chains in conjunction with the OCFE method, the achieved computational gain is balanced with the loss of accuracy on the MWD calculations.