Population Balance Model (PBM)

MR-CFD experts are ready for PBM analysis, consulting, training, and CFD simulation.

Population Balance Model (PBM)

The Population Balance Model (PBM) is one of the concepts in the field of collective phenomenology. This concept is used when components or particles move through state space or during birth and death processes, the current particles are destroyed, and then new particles are produced.

The main application of this model is related to currents that include a secondary phase with a specific size distribution. Particle size distribution, including solid particles, bubbles, droplets, etc., can evolve due to transfer within the base fluid or due to the occurrence of chemical reactions in a multiphase system. Evolutionary processes include nucleation, growth, Aggregation, and Breakage.


Consider, for example, a precipitation process. This process is related to forming solid particles from a liquid solution and consists of several processes, including nucleation, growth, Aggregation, failure, etc. These processes result in an increase or decrease in the number of particles of a specific size. So it can be said that PBM means balance in the number of particles with a particular state (for example, particle size in the example).

PBM Application

The issue of particle PBM is used in various topics, including crystallization processes, precipitative reactions, bubble columns, sprays, polymerization with fluidized bed, grain Granulation, liquid-liquid emulsion, aerosol flows, leaching metallurgy, liquid-liquid extraction, gas-liquid dispersion -liquid dispersions), liquid-liquid reactions, biology, etc.

Population Balance Equation Application

A PBM equation is needed to describe population or particle density changes in multiphase flows involving a given size distribution. This balance equation is commonly known as the PBM equation. In fact, in multiphase currents with particle distribution and the equilibrium equations of momentum, mass, and energy, the population balance equation must also be solved.

Population Balance Equations are defined in several branches of modern science (mainly chemical engineering). These equations are used to describe the gradual evolution of particles and the motion of a particle population. In other words, these equations define how individual particle populations develop over time in a series of distinct thermophysical properties. These equations include a set of integral-partial differential equations (PDE), each of which presents the field behavior of a population of particles based on the behavior of each independent particle in the local state.

To use this concept in modeling, it is necessary to introduce the “number density” function and the “state vector” to calculate the population or particle density. Different particles can be distinguished from the particle population. Their behavior can be described using these concepts and with the help of particle properties (such as particle size, structural composition, porosity, etc.).

Number Density

The number density is considered as an intensity quantity. It is equal to the number of countable components (such as the number of particles, molecules, atoms, photons, etc.) per unit volume. In other words, the number density is considered a criterion for describing the degree of concentration of the desired particles. Of course, the number density can be volumetric (three-dimensional), surface (two-dimensional), and linear (one-dimensional).

The number density per unit volume is expressed as the following relation:


N equals the number of particles, and V equals the volume studied in the above relation. Also, the unit of number density is equal to 1 / m3. It should be noted that the value of N should be so significant that if the resulting value of n is rounded to the nearest integer, no problem will appear. The value of V should be so small that the value of n does not profoundly depend on the size or shape of volume V.

When the number density is defined as a function of spatial coordinates, the expression N is equal to the total number of particles in the whole volume V:


State Vector

The state vector of a particle, by a series of external coordinates (for example, vector x), expresses the particle’s location and internal coordinates (with the symbol φ) to express items such as particle size, particle composition, and particle temperature, etc. In the present discussion, the internal coordinates include only the particle size. The state vector of each particle with the symbol s (x, φ) depends on the local occupied volume around the particle (according to the spatial coordinates of the particle) and the volume of the particle itself (proportional to the particle size).

According to these internal and external coordinates, the function of the number density is expressed as follows:


So the total number of particles throughout the system is equal to:


By considering the amount of local volume around the desired spatial coordinates, the average number density value in the desired spatial space is obtained:


According to the above relation, by dividing the number of particles by the amount of local volume around the particle, the average number density value is obtained. Now, assuming a certain amount of volume for each particle (according to the particle size) and multiplying this particle volume by the sides of the above relation, the expression related to the volume fraction of all particles is obtained:


Population Balance Equations

Assuming that φ is equivalent to a certain amount of volume (V), the Transport Equation is expressed as follows:


It should be noted that in the above equations, the time term is also considered.

The above equations include the terms growth rate, birth due to Aggregation, death due to Aggregation, birth due to Breakage, and death due to Breakage.



Also, the boundary conditions and initial conditions are stated as follows:


Whereas in the above relation, the expression  indicates the nucleation rate in particles with number/m3.s.

Particle Growth

Particle growth rate based on particle volume and particle diameter growth rate (same as length dimension) is defined as follows:


Also, the volume of a single particle (equivalent to V) is defined as follows, and as a result, the relationship between Gv and GL is written as follows:



The phenomena of birth and death can occur as a result of the process of Breakage. Examples of the breakage process are the Breakage of the crystal during the crystallization process and the Breakage of the bubble due to the turbulence of the liquid inside a liquid column.

The breakage rate term is expressed as follows:


In the above relation, g (V’) equals the Breakage frequency (meaning a fraction of particles with a volume of V’ breaking per unit time per m3/s) and β (V│V’) equals the probability density function (PDF) of particles which are breaking from a volume V’ to a particle volume equal to V.

The birth rate of particles with a volume V due to failure is expressed as follows:


Thus, the number equal to particles of volume V’, breaks per unit time to produce the number of  particles. P is also equal to the number of produced parent particles (for example, two particles due to Breakage).

The particle death rate with volume V due to Breakage is expressed as follows:


It should be noted that in statistics and probability, the probability density function (PDF) of a hypothetical dependent variable (x) in a given interval (x1, x2) is equivalent to the integral value of that function in the same specified interval of the dependent variable. So the value of the integral represents the probability density function:


In this case, β (V│V’), as a probability density function (PDF), represents the function of the fragmentation distribution of particles or the size distribution of productive particles; That is, the integral of the density function of the number of particles in a given range of particle size is equivalent to the probability function of the desired particle size.

Breakage in ANSYS Fluent

To define the breakage process in ANSYS Fluent software, defining the frequency (in units of 1/s) and the PDF is necessary. To define a frequency, you can specify a constant value or define a UDF based on available equations. The software also has several models for estimating the frequency, including Luo, Lehr, and Ghadiri models. However, the software also uses various models for the PDF, including the parabolic, Laakkonen, and generalized models.

Luo and Lehr’s models integrate a term that simultaneously considers the breakage frequency and the PDF of the breakage particles. This means that by selecting one of the two methods of determining the breakage frequency, there is no need to specify the PDF, and it is automatically defined with the breakage frequency function. The equation for these two methods of calculating the breakage frequency is as follows:


So ξ equals λ/d, which means the ratio of the size of the striking vortices (λ) to particles with diameter d. Some of the parameters of the above equation are presented in the table below:

m b n K (1/m3.s) model
-11/3 12(ƒ2/3+(1-ƒ)2/3-1)σρ-1ε-2/3d-5/3 11/3 0.9238ε1/3d-2/3α Luo
-2/3 2σρ-1ε-2/3d-5/3ƒ-1/3 13/3 1.19ε-1/3d-7/3σρ-1ƒ-1/3 Lehr


The Ghadiri model, on the other hand, only models the frequency of Breakage solid particles, and the PDF model must be determined independently to determine the distribution of productive particles. The equation for this method is to calculate the breakage frequency related to the material properties and the particle impact state and is expressed as follows:


So ρs represents the particle density, E represents the elastic modulus, T represents the interfacial energy, v represents the impact velocity, and L represents the particle diameter before Breakage.

So for this model, a breakage PDF must be used. The parabolic model of the PDF function contains information about the probability of parts formed due to a breakage phenomenon, as well as the number of particles and the possible size distribution due to Breakage. This parabolic model defines the PDF related to Breakage as follows:


Depending on the value of the shape factor with the symbol C, different behaviors will be observed in the shape of the particle Breakage distribution function. For example, suppose C is equal to 2. In that case, the particle Breakage has a uniform distribution. If C is between zero and 2, a concave parabola (parts of unequal size) is obtained, and if C is between 2 and 3, A convex parabola (parts of unequal size) is obtained.

The generalized PDF of the Breakage function simulates the components resulting from multiple Breakages (Breakage process in more than two steps) and the form of distribution of productive particles (including uniform, uniform size, erosion, power law, parabolic, and dual beta). The generalized model of the PDF is presented as follows:


In the above relation, wi represents the weight factor, pi represents the average number of productive particles, and qi and ri represent the power representation.


Similar to Breakage, the phenomena of birth and death can also occur due to Aggregation. Examples of the aggregation process are the composition of particles during the crystallization process and the bubble fusion within the liquid column reactors.

The expression of the aggregation rate is expressed as below a(V, V’). The aggregation rate term has a unit of m3/s. It is defined as the product of two quantities, including the collision frequency between particles with volume V and particles with volume V’, and the efficiency of Aggregation (meaning the probability of fusion of particles of volume V with of V’).

The particle birth rate with a volume V due to the aggregation process is defined as follows:


While particles with a volume of V-V’combine with particles with a volume of V’ to form particles with a volume of V. A factor of 1/2 is also used to prevent double-counting of each stroke.

The particle death rate with a volume V due to the aggregation process is defined as follows:


It should be noted that the Breakage and Aggregation rates depend on the nature of the physical application of the problem. For example, in gas-liquid dispersion, the Breakage and Aggregation rates are considered local gas-liquid turbulence loss functions.

Aggregation in ANSYS Fluent

A constant value can be set to define the aggregation process in Ansys Fluent software, or a UDF can be defined based on the available equations. The software also has several models for estimating aggregation rates, including Luo, free molecular, and turbulent models.

According to the Luo model, the rate of Aggregation is defined as the rate of particle volume formation due to double collisions of particles with volumes Vi and Vj:


In the above relation, is equivalent to the collision frequency, and  is equivalent to the probability that the collision will lead to confusion. Aggregation Frequency is defined as the following relation:


Where uij represents the velocity of two particles colliding with the diameters di and dj and the densities of si and sj.

According to the free molecular aggregation model, real particles aggregate and break at a frequency determined by complex dependencies on the internal coordinates of the particles. Typically, tiny particles (up to one micrometer) aggregate due to collisions caused by Brownian motion. For example, the collisions frequency is size-dependent and generally follows the following relation:


So that kB is equal to the Boltzmann constant and T is equal to the absolute temperature and µ is equal to the viscosity of the suspended fluid.

The turbulent model is such that mechanical energy is transferred to the fluid during the mixing process, and this energy causes turbulence within the fluid. Turbulence creates vortices that lead to energy waste. Energy is transferred from the largest vortices to the smallest vortices, within which energy is lost through viscous interactions. The size of the smallest vortices is equivalent to the Kolmogorov microscale, which is expressed as a function of kinematic viscosity and turbulence energy dissipation rate:


In turbulent flow fields, the aggregation process can take two mechanisms: the viscous and inertial mechanisms. The viscosity mechanism is applied when the particles are smaller than the Kolmogorov microscale. The mechanism of inertia, however, applies when the particles are larger than the Kolmogorov microscale. In this case, it is assumed that the particles are independent of the velocity values.

For the viscous mechanism, particle collisions are affected by local shear within the vortex. Thus, the aggregation rate is expressed as follows:


For the inertia mechanism, the particles are larger than the smallest vortex, and hence, the particles are pulled by fluctuations in velocity within the flow field. In this case, the aggregation rate is expressed as follows:


Birth Due to Nucleation

Due to the application of the study environment, the spontaneous nucleation process of particles can occur due to the transfer of molecules from the initial phase. For example, in forming crystals from solution, the first step is related to the phenomenon of fuzzy separation or the birth of new crystals. In welding applications, the formation of elementary vapor bubbles is considered a nucleation process called the nuclear boiling process. The nucleation rate is defined by the boundary condition shown in the following statement:


Solution Methods

ANSYS Fluent software can solve population balance equations using several different methods. These solving population balance equations include the Discrete Method, the standard method of moments (SMM), and the quadrature method of moments (QMOM). For each of these methods, the execution of the software is limited to a single internal coordinate associated with the particle size. Also, each of these methods has advantages and disadvantages.

Discrete Method

The discrete method is based on a continuous particle size distribution (PSD) representation, which acts in a series of size classifications (bins). The advantage of this method is its high power of the numerical solution, and the particle size distribution can be obtained directly. However, the disadvantage of this method is that bins must be obtained by inductive reasoning; This means that by solving the model in a few steps, a suitable number of required categories must be achieved. Another disadvantage of this method is that sometimes many bins may be required, in which case the software solution process is slow.

The following image shows a particle size distribution based on the discrete method:


According to this method in ANSYS Fluent software, the particle balance equation (PBE) is defined in terms of particle size fraction i volume fraction terms as follows:


So ρs represents the density of the secondary phase and αi as the volume fraction of the particle size i. Thus:


And also


In the above equation, there are the terms of nucleation rate, growth rate, and birth and death rates due to Aggregation or Breakage. For example, the nucleation rate is shown as  for the smallest volume fraction (V0) in the discrete equation. The sign 0i also indicates that this particle term () in the above equation appears only if the smallest particle size is desired.

The growth rate in the above equation is expressed as follows:


Then the birth and death rates due to the phenomena of Aggregation and Breakage are expressed in the following terms:


In the discrete method, the breakage formulation is Hagesather and Ramakrishna. According to the Hagesather method, the breakage sources are distributed in sequential size categories, maintaining the mass and numerical density. For the case that the ratio between the sizes of consecutive categories can be expressed in the form (s=1,2,…)2s, the source in the i category (i = 1,… .., n) can be expressed as follows:


However, the Ramakrishna method is considered as a more accurate formulation, and according to this method, the Breakage rate is expressed as follows:


The Ramakrishna formulation can be slow due to the high number of integral points required. However, integrals for simple forms of β can be relatively easy to perform.

It should be noted that the Hagesather formulation requires fewer integral points, and the weakness of this method compared to the Ramakrishna formulation in terms of accuracy can be corrected by choosing the right size for batch sizes.

Standard Method of Moments (SMM)

The standard moment method (SMM) is an alternative method for solving the population balance equation (PBE). One of the advantages of this method is that it reduces the problem in terms of dimensions. Secondly, this method is relatively simple to solve the transport equations with a lower degree of momentum. However, this method also has disadvantages; Thus, firstly, the Aggregation phenomenon must be constant, and growth must be independent of particle size. Secondly, it is not possible to model the phenomenon of Breakage.

The solution of the SSM method is based on the momentum of the population balance equation (PBE) concerning the internal coordinates (in this case, the particle size is L as internal coordinates). So the moment k is defined as follows:


Assuming a constant value for particle growth, the transport equation can be written as follows:


In the above expressions, n is equivalent to a certain number of moments and is equivalent to the nucleation rate. The growth term is also written as follows:


Quadrature Method of Moments (QMOM)

The quadrature method of moments (QMOM) models the gradual motion and deformation of sprays and coagulation problems. This method requires a relatively small number of scalar equations to follow population moments with minor errors. This method is an excellent alternative to the discrete method when aggregate values ​​are desired instead of an exact particle size distribution. Features of this method include the need for fewer variables (usually only six or eight moments) and performing a dynamic calculation for size categories. In contrast, the disadvantage of this method is that the number of horizontal distances (distance of points from the vertical axis or abscissas) may not be sufficient to describe the population size distribution (PSD) and that the solution of the product-difference algorithm may be time-consuming.

The quadrature approximation is based on the sequence of orthogonal polynomials with s(L) (e.g., particle size distribution). Suppose the horizontal distance of the quadrature approximation is the same as the polynomial nodes of degree n. In that case, the quadrature approximation written in the following statement is exact, provided that f(L) is a polynomial of degree n or less.


A direct way to calculate the quadrature approximation is by defining it through moments:


The quadrature approximation of degree n is defined by n to the weight of wi and n to the horizontal distance Li and can be denoted by 2n to its first moment, including m0,…..,m2n-1 recursive relations for polynomials as Mk moments are calculated. Once this relation is written in the form of a matrix, the polynomial roots associated with specific values ​​of the Jacobi matrix can be easily represented. This process is known as the product-difference algorithm. Once the horizontal weights and distances have been determined, the source terms due to Aggregation and Breakage can be calculated, and hence, the transport equations for the moments can be solved. Therefore, the terms of birth and death in the above equation are written in the following forms:


In general, in theory, there is no limit to the expression of breakage and aggregation rates when using the quadrature method of moment (QMOM).

Also, the nucleation rate is similarly defined by the standard moment method (SMM). However, the growth rate for the QMOM is defined as follows:


It should be noted that the QMOM is possible for the growth rate depending on the size.


MR-CFD experts are ready to fulfill every Computational Fluid Dynamic (CFD) needs. Our service includes both industrial and academic purposes considering a wide range of CFD problems. MR-CFD services in three main categories of ANSYS Fluent Consultation, ANSYS Fluent Training, and ANSYS Fluent Project Simulation. MR-CFD company has gathered experts from various engineering fields to ensure the quality of CFD services. Your CFD project would be done in the shortest time, with the highest quality and appropriate cost.

Service Process

Here is MR-CFD Service Process Steps


Contact Us via [email protected] or Call on WhatsApp to Share the Project Description.


Project Order Will Be Investigated by MR-CFD Experts.


An Official Contract Including Service Price, Time Span, and Terms of Condition Will Be Set as a Service Agreement.


All Simulation Files, Results, Technical Report, and a Free Comprehensive Training Movie Will Be Sent to the Client as the Contract is Done.

What are you waiting for? Submit your project consulting for free!

Free Consulting!

Yes, it is right that we can do everything you want the only thing that you should do is to fill the form at the Free Consulting Page and send us your project’s details then our support team will be in touch with you as soon as possible and they will give you free consulting for your project.

Free Consulting
Back To Top

Refund Reason

you tube
Call On WhatsApp