Numerical solution of population balance equations for nucleation, growth and aggregation processes

This article focuses on the derivation of numerical schemes for solving population balance models (PBMs) with simultaneous nucleation, growth and aggregation processes. Two numerical methods are proposed for this purpose. The first method combines a method of characteristics (MOC) for growth process with a finite volume scheme (FVS) for aggregation process. For handling nucleation terms, a cell of nuclei size is added at a given time level. The second method purely uses a semi-discrete finite volume scheme for nucleation, growth and aggregation of particles. Note that both schemes use the same finite volume scheme for aggregation process. On one hand, the method of characteristics offers a technique which is in general a powerful tool for solving linear growth processes, has the capability to overcome numerical diffusion and dispersion, is computationally efficient, as well as give highly resolved solutions. On the other hand, the finite volume schemes which were derived for a general system in divergence form, are applicable to any grid to control resolution, and are also computationally not expensive. In the first method a combination of finite volume scheme and the method of characteristics gives a highly accurate and efficient scheme for simultaneous nucleation, growth and aggregation processes. The second method demonstrates the applicability, generality, robustness and efficiency of high-resolution schemes. The proposed techniques are tested for pure growth, simultaneous growth and aggregation, nucleation and growth, as well as simultaneous nucleation, growth and aggregation processes. The numerical results of both schemes are compared with each other and are also validated against available analytical solutions. The numerical results of the schemes are in good agreement with the analytical solutions.     

On the solution of population balance equations (PBE) with accurate front tracking methods in practical crystallization processes

The PBE (population balance equation) containing birth, growth, agglomeration and breakage kinetics is described by a conservation law with a moving source term. For the solution of the PBE, we compare two accurate front tracking methods such as a modified method of characteristics (MOC) and a finite difference method with the weighted essentially non-oscillatory (WENO) scheme. Both methods are applied to a potassium sulfate crystallization problem (K₂SO₄-H₂O system) with a discontinuous initial condition. Parameters of agglomeration and breakage kinetics are estimated on the basis of the experimental data of the K₂SO₄-H₂O system.Owing to moving axis along a crystal growth rate (i.e. elimination of the growth term), the modified MOC is able to provide a highly accurate solution even at discontinuous points without numerical diffusion error. However, in the case of stiff nucleation which can commonly appear in practical crystallization processes, it is necessary to adaptively determine time levels to add a mesh of the nuclei size. For solving PBEs involving agglomeration and breakage terms, the MOC can take more long computational time than the spatial discretization methods like the WENO scheme. It is pointed out that the MOC is not available to solve more than two coupled PBEs in general.WENO schemes for spatial discretization are firstly addressed in this study for the dynamic simulation of batch crystallization processes. The WENO schemes show improvements of accuracy and stability over conventional discretization methods (e.g., backward, central or common upwinding schemes). However the WENO schemes on fixed meshes show, to some extent, the numerical diffusion error near discontinuities or steep moving fronts like other finite difference methods. Hence, they require spatially-adaptive mesh techniques in order to track more accurately the moving fronts. Even though the WENO schemes are less accurate than the MOC, they are of practical use for solving complex PBEs owing to a short computational time and little limitation to use.     

Comparison of numerical methods for the simulation of dispersed phase systems

This contribution deals with a new numerical method for an accurate and efficient simulation of particulate processes. As an example for dispersed phase systems a detailed model for crystallization processes is considered. After the model derivation, which is based solely on physical principles, different techniques for the numerical simulation of the mathematical model are discussed. State-of-the-art finite volume schemes based on the ‘method of lines’ approach are compared to the recently published ‘Space-time conservation element and solution element method’. The presented simulation results show a strong dependence on the chosen numerical method. Guidelines for a proper selection of numerical methods for the treatment of population balance based models are given.     

Model based robust control approach for batch crystallization product design

The paper presents a novel control approach for crystallization processes, which can be used for designing the shape of the crystal size distribution to robustly achieve desired product properties. The approach is based on a robust optimal control scheme, which takes parametric uncertainties into account to provide decreased batch-to-batch variability of the shape of the crystal size distribution. Both open-loop and closed-loop robust control schemes are evaluated. The open-loop approach is based on a robust end-point nonlinear model predictive control (NMPC) scheme which is implemented in a hierarchical structure. On the lower level a supersaturation control approach is used that drives the system in the phase diagram according to a concentration versus temperature trajectory. On the higher level a robust model-based optimization algorithm adapts the setpoint of the supersaturation controller to counteract the effects of changing operating conditions. The process is modelled using the population balance equation (PBE), which is solved using a novel efficient approach that combines the quadrature method of moment (QMOM) and method of characteristics (MOC). The proposed robust model based control approach is corroborated for the case of various desired shapes of the target distribution.     

On the solution of population balances for nucleation, growth, aggregation and breakage processes

This work is concerned with the modeling and simulation of population balance equations (PBEs) for combined particulate processes. In this study a PBE with simultaneous nucleation, growth, aggregation and breakage processes is considered. In order to apply the finite volume schemes (FVS) a reformulation of the original PBE is introduced. This reformulation not only help us to treat the aggregation and breakage processes in a manner similar to the growth process in the FVS but also in deriving a stable numerical scheme. Two numerical methods are proposed for the numerical approximation of the resulting reformulated PBE. The first method combines a method of characteristics (MOC) for growth process with an FVS for aggregation and breakage processes. The second method purely uses a semidiscrete FVS for all processes. Both schemes use the same FVS for aggregation and breakage processes. The numerical results of the schemes are compared with each other and with the available analytical solutions. The numerical results were found to be in good agreement with analytical solutions.     

An efficient numerical technique for solving one-dimensional batch crystallization models with size-dependent growth rates

In this work, an efficient numerical method is introduced for solving one-dimensional batch crystallization models with size-dependent growth rates. The proposed method consist of two parts. In the first part, a coupled system of ordinary differential equations (ODEs) for the moments and the solute concentration is numerically solved to obtain their discrete values in the time domain of interest. These discrete values are also used to get growth and nucleation rates in the same time domain. To overcome the issue of closure, a Gaussian quadrature method based on orthogonal polynomials is employed for approximating integrals appearing in the ODE system. In the second part, the discrete growth and nucleation rates along with the initial crystal size distribution (CSD) are used to construct the final CSD. The expression for CSD is obtained by applying the method of characteristics and Duhamel's principle on the given population balance model (PBM). The proposed method is efficient, accurate, and easy to implement in the computer. Several numerical test problems of batch crystallization processes are considered. For a validation, the results of the proposed technique are compared with those obtained using a high resolution finite volume scheme.     

A control oriented study on the numerical solution of the population balance equation for crystallization processes

The population balance equation provides a well established mathematical framework for dynamic modeling of numerous particulate processes. Numerical solution of the population balance equation is often complicated due to the occurrence of steep moving fronts and/or sharp discontinuities. This study aims to give a comprehensive analysis of the most widely used population balance solution methods, namely the method of characteristics, the finite volume methods and the finite element methods, in terms of the performance requirements essential for on-line control applications. The numerical techniques are used to solve the dynamic population balance equation of various test problems as well as industrial crystallization processes undergoing simultaneous nucleation and growth. The time-varying supersaturation profiles in the latter real-life case studies provide more realistic scenarios to identify the advantages and pitfalls of a particular numerical technique.The simulation results demonstrate that the method of characteristics gives the most accurate numerical predictions, whereas high computational burden limits its use for complex real crystallization processes. It is shown that the high order finite volume methods in combination with flux limiting functions are well capable of capturing sharp discontinuities and steep moving fronts at a reasonable computational cost, which facilitates their use for on-line control applications. The finite element methods, namely the orthogonal collocation and the Galerkin's techniques, on the other hand may severely suffer from numerical problems. This shortcoming, in addition to their complex implementation and low computational efficiency, makes the finite element methods less appealing for the intended application.     

FE coupled to SPH numerical model for the simulation of high-velocity impact on ceramic based ballistic shields

Predictive models are an important tool in the design and optimization of ballistic shields. Indeed, several authors in the literature have developed numerical models for simulating high-velocity impact on ceramic-based ballistic shields which are based on the finite element method. Element erosion is usually implemented in finite element models simulating impact to remove excessively distorted elements but, it leads to energy loss, which in turns may lead to the production of incorrect results. Due to the absence of a fixed mesh, the smoothed particle hydrodynamics method is well suited for large deformation problems, overcoming the limitations of the finite element method. On the other hand, the smoothed particle hydrodynamics method is computationally more expensive than the finite element method. Thus, a numerical model combining the lower computational cost of finite elements and the capability of smoothed particle hydrodynamics of dealing with crack formation and fracturing would be an interesting solution for the simulation of high-velocity impact on ceramics. The aim of this work is therefore to develop a finite element coupled to smoothed particle hydrodynamics numerical model for the simulation of high-velocity impact on ceramic-based ballistic shields. High-velocity impact tests were performed on Al₂O₃ tiles and the experimental results were used for the calibration of the numerical models; furthermore, high-velocity impact test were performed on multilayer targets with Al₂O₃ front layer and AA6061-T6 backing layer for the validation of the numerical models. This study proved that this approach is more appropriate for the simulation of the response of ceramic materials rather the common finite element model.     

A Laplace transformation based technique for reconstructing crystal size distributions regarding size independent growth

This article introduces a technique for reconstructing crystal size distributions (CSDs) described by well-established batch crystallization models. The method requires the knowledge of the initial CSD which can also be used to calculate the initial moments and initial liquid mass. The solution of the reduced four-moment system of ordinary differential equations (ODEs) coupled with an algebraic equation for the mass gives us moments and mass at the discrete points of the given computational time domain. This information can be used to get the discrete values of size independent growth and nucleation rates. The discrete values of growth and nucleation rates along with the initial distribution are sufficient to reconstruct the final CSD. In the derivation of current technique the Laplace transformation of the population balance equation (PBE) plays an important role. The proposed technique has dual purposes. Firstly, it can be used as a numerical technique to solve the given population balance model (PBM) for batch crystallization. Secondly, it can be used to reconstruct the final CSD from the initial one and also vice versa. The method is very efficient, accurate and easy to implement. Several numerical test problems of batch crystallization processes are considered here. For validation, the results of the proposed technique are compared with those from the high resolution finite volume scheme which solves the given PBM directly.     

Part IV: Dynamic evolution of the particle size distribution in particulate processes. A comparative study between Monte Carlo and the generalized method of moments

The present work provides a comparative study on the numerical solution of the dynamic population balance equation (PBE) for batch particulate processes undergoing simultaneous particle aggregation, growth and nucleation. The general PBE was numerically solved using three different techniques namely, the Galerkin on finite elements method (GFEM), the generalized method of moments (GMOM) and the stochastic Monte Carlo (MC) method. Numerical simulations were carried out over a wide range of variation of particle aggregation and growth rate models. The performance of the selected techniques was assessed in terms of their numerical accuracy and computational requirements. The numerical results revealed that, in general, the GFEM provides more accurate predictions of the particle size distribution (PSD) than the other two methods, however, at the expense of more computational effort and time. On the other hand, the GMOM yields very accurate predictions of selected moments of the distribution and has minimal computational requirements. However, its main disadvantage is related to its inherent difficulty in reconstructing the original distribution using a finite set of calculated moments. Finally, stochastic MC simulations can provide very accurate predictions of both PSD and its corresponding moments while the MC computational requirements are, in general, lower than those required for the GFEM.     

Part I: Dynamic evolution of the particle size distribution in particulate processes undergoing combined particle growth and aggregation

The present study provides a comprehensive investigation on the numerical problems arising in the solution of dynamic population balance equations (PBEs) for particulate processes undergoing simultaneous particle growth and aggregation. The general PBE was numerically solved in both the continuous and its equivalent discrete form using the orthogonal collocation on finite elements (OCFE) and the discretized PBE method (DPBE), respectively. A detailed investigation on the effect of different particle growth rate functions on the calculated PSD was carried out over a wide range of variation of dimensionless aggregation and growth times. The performance (i.e., accuracy and stability) of the employed numerical methods was assessed by a direct comparison of predicted PSDs or/and their respective moments to available analytical solutions. It was found that the OCFE method was in general more accurate than the discretized PBE method but was susceptible to numerical instabilities. On the other hand, for growth dominated systems, the discretized PBE method was very robust but suffered from poor accuracy. For both methods, discretization of the volume domain was found to affect significantly the performance of the numerical solution. The optimal discretization of the volume domain was closely related with the satisfactory resolution of the time-varying PSD. Finally, it was shown that, in specific cases, further improvement of the numerical results could be obtained with the addition of an artificial diffusion term or the use of a moment-weighting method to correct the calculated PSD.     

Integration and dynamic inversion of population balance equations with size-dependent growth rate

Two novel solution schemes for integration and dynamic inversion of a class of population balance equations with size-dependent growth rate are contributed in this article. The proposed methods are developed for population balance systems that, in addition to an advective and a birth rate term, include an external variable which may be fixed (referring to the integration problem), or is to be computed for some prespecified evolution of the population system (referring to the dynamic inversion problem). A unique diffeomorphism for the independent time and internal (property) coordinates is introduced, which transforms the original nonlinear partial differential equation into a linear one with straight lines as equation characteristics. The evolution of the density function in the time and internal property coordinates is then computed by transporting the initial and boundary density functions in the transformed domain. For the integration of the temporal behavior of the boundary density function, a generalization of the standard method of moments is introduced, resulting in a closed integro-differential structure driven by convolution and correlation integrals. While the correlation integrals refer to the given initial density function and, hence, can be a priori computed, the convolution integrals involve the boundary density function and have to be integrated a posteriori online. The solution of the dynamic inversion problem, on the other hand, turns out to become simpler, as it converts to an algebraic equation after pre-computation of the correlation/convolution integrals for the given initial and boundary density function. In a next step, we introduce the concept of internal or eigenmoments which is useful for the representation of the original physical moments in terms of an infinite series. The dynamics of eigenmoments exhibits a closed ODE structure, which we refer to as the internal model. From the perspective of the systems theory this turns out to be an infinite dimensional flat system. Hence, in addition to providing a highly simple structure – basically, a chain of integrators – for the integration of the moments and the density function, the internal model allows for a direct solution of the dynamic inversion problem due to its flatness property. However, the ease and elegance of the method, in general, come at the price of an approximation of the infinite dimensional problem by a finite one. The usability of both proposed solution methods is illustrated on a batch crystallization process with size-dependent growth rate kinetics. The proposed methods are compared in terms of efficiency and accuracy with a state-of-the-art high-resolution finite volume scheme using a numerical example.     

熔体充模过程动态模拟及流场分析

将Level Set／Ghost方法应用于聚合物成型研究，实现了非等温情况下注射成型聚合物熔体充模阶段的动态模拟；得到了正确的流线分布和不同时刻的温度、压力等值线分布。对Level Set／Ghost方程的求解，空间方向采用高分辨率、稳定且无振荡的5WENO（the fifth—order weighted essentially non—oscillatory）格式进行离散，时间方向采用稳定的TVD-RK（total variation diminishing Runge Kutta）方法进行离散。物理量控制方程采用一般的有限差分格式进行数值求解。结果表明，Level Set／Ghost方法可以准确追踪非等温聚合物熔体前沿界面的位置，并能精确描述前沿界面的形状，同时可以实现动态流场物理量的准确模拟。     

Simulation of rapid pressure swing adsorption and reaction processes

A general model for non-isothermal adsorption and reaction in a rapid pressure swing process is described. Several numerical discretisation methods for the solution of the model are compared. These include the methods of orthogonal collocation, orthogonal collocation on finite elements, double orthogonal collocation on finite elements, and cells-in-series. Computationally, orthogonal collocation on finite elements is found to be the most efficient of these. The model is applied to air separation for oxygen production. Calculations confirm the formation of a concentration shock when an adsorbent bed is pressurised with air. The form and propagation of the shock over short times is found to be in excellent agreement with the exact similarity transformation solutions derived for an infinitely long bed. For air separation, novel experimental measurements, showing an optimum particle size for maximum product oxygen purity, are accurately described by the model. Calculations indicate that a poor separation results from ineffective pressure swing for beds containing very small particles, and from intraparticle diffusional limitations for beds containing very large particles. For adsorption coupled with reaction, finite rate and reversible reactions are considered. These include both competitive and non-competitive reaction schemes. For the test case of a dilute reaction A &.rlhar2; B + 3C, with B the only adsorbing species, bed pressurisation calculations are found to be in excellent agreement with the solutions obtained by the method of characteristics.     

Traitement par une méthode d'éléments finis de modèles de colonnes de rectification discontinue à garnissage

Mathematicals models representing several operating modes of a packed batch distillation column are detailed. In the first part, the authors consider steady-state operation with total reflux, where the axial and radial dispersion phenomena are superimposed on the plug flow of the liquid and vapour phases. The equations are solved by means of a finite element method based on the Galerkin criterion, and also by using a finite difference procedure. The numerical study of the steady-state operating mode, including both axial and radial dispersion phenomena, allows us to note that the last term can be, in most cases, neglected. On the other hand, the dynamic model analysis has shown that the axial dispersion phenomenon is always indispensable for accurately representing the column operation.     

Caractérisation viscoélastique du comportement d'une membrane thermoplastique et modélisation numérique de thermoformage

In this work, we are interested, on the one hand in the characterization of circular polymeric ABS membrane under biaxial deformation using the bubble inflation technique, on the other hand in modelling and numerical simulation of the thermoforming of ABS materials using the dynamic finite element method. The viscoelastic behaviour of the Lodge model is considered. First, the governing equations for the inflation of a flat circular membrane are solved using a variable‐step‐size‐finite difference method and a modified Levenberg‐Marquardt algorithm to minimize the difference between the calculated and measured inflation pressure. This will determine the material constants embedded within the model used. For dynamic finite elements method, we consider a nonlinear load in air flow which obeys the Redlich‐Kwong equation of state of the real gases. For numerical simulation, the lagrangian formulation together with the assumption of the membrane theory is used. Moreover, the influence of the viscoelastic model on the thickness and on the stress distribution in the thermoforming sheet are analysed for ABS material.     

Parallel and high resolution numerical solution of concentration and temperature distributions in fluidized beds

In this article higher order in time numerical schemes with efficient time stepping for the solution of concentration and temperature distributions in fluidized beds using parallel computers are presented. The mathematical model equations consist of strongly coupled and semi linear convection-diffusion-reaction equations. Invariant regions for the model are derived to check the solution bounds. The numerical discretization for the space using the finite element method is presented and the numerical treatment is enhanced by using adaptive and higher order linearly implicit Runge–Kutta methods for the time discretization. For different time stepping methods and different spatial grid sizes numerical results are obtained and compared. The methods used show a clear improvement for the problem under consideration compared to previously presented results (Nagaiah, Warnecke, Heinrich, & Peglow, 2008). Additionally, the higher order time stepping methods yield a good parallel efficiency, paving the way for the efficient study of more complex phenomena.     

Mathematical simulation of wet spinning coagulation process: Dynamic modeling and numerical results

A dynamic model of polymer wet spinning coagulation process is proposed in this article. The model is based on the double diffusion phenomenon, phase separation process, continuity balance, and momentum balance of the entire coagulation process. The uniqueness of the model lies in its dynamic feature. The model can simulate the system's dynamic response to variations in system inputs/parameters. Steady‐state system solutions can also be produced as the long‐time solutions of the dynamic model; a settling time can be observed at the same time. This paper employs a computationally efficient method of lines numerical algorithm for solving the dynamic model. A simulation experiment on a selected non‐solvent‐solvent‐polymer ternary system is carried out to verify the model as well as the numerical method. The dynamic simulation results are analyzed and discussed. At the end of the article, h‐refinement and p‐refinement are used to confirm the spatial convergence of the numerical solutions. © 2016 American Institute of Chemical Engineers AIChE J, 62: 3432–3440, 2016     

High‐order simulation of polymorphic crystallization using weighted essentially nonoscillatory methods

Most pharmaceutical manufacturing processes include a series of crystallization processes to increase purity with the last crystallization used to produce crystals of desired size, shape, and crystal form. The fact that different crystal forms (known as polymorphs) can have vastly different characteristics has motivated efforts to understand, simulate, and control polymorphic crystallization processes. This article proposes the use of weighted essentially nonoscillatory (WENO) methods for the numerical simulation of population balance models (PBMs) for crystallization processes, which provide much higher order accuracy than previously considered methods for simulating PBMs, and also excellent accuracy for sharp or discontinuous distributions. Three different WENO methods are shown to provide substantial reductions in numerical diffusion or dispersion compared with the other finite difference and finite volume methods described in the literature for solving PBMs, in an application to the polymorphic crystallization of L ‐glutamic acid. © 2008 American Institute of Chemical Engineers AIChE J, 2009