Analytical Approach For Solving Population Balances: A Homotopy Perturbation Method
Gurmeet Kaur Randhir Singh Mehakpreet Singh Jitendra Kumar Themis Matsoukas
math.NA




Analytical Approach For Solving Population Balances: A Homotopy Perturbation Method



Gurmeet Kaur^a, Randhir Singh^b, Mehakpreet Singh^c, Jitendra Kumar^a, Themis Matsoukas^d




^aDepartment of Mathematics, Indian Institute of Technology Kharagpur, Kharagpur, India

^bDepartment of Mathematics, Birla Institute of Technology  Mesra, Ranchi-835215, India

^cBIOMATH, Department of Mathematical Modelling, Statistics and Bioinformatics,

Faculty of Bioscience Engineering, Ghent University, 9000 Ghent, Belgium 

^dDepartment of Chemical Engineering, Pennsylvania State University, 158 Fenske Laboratory,
University Park, PA 16802-4400, USA







In the present work, a new approach is proposed for finding the analytical solution of population balances. This approach is relying on idea of Homotopy Perturbation Method (HPM). The HPM solves both linear and nonlinear initial and boundary value problems without nonphysical restrictive assumptions such as linearization and discretization. It gives the solution in the form of series with easily computable solution components. The outcome of this study reveals that the proposed method can avoid numerical stability problems which often characterize in general numerical techniques related to this area. Several examples including Austin's kernel available in literature are examined to demonstrate the accuracy and applicability of the proposed method.


Keyword: Particles; Population Balance Equation; Homotopy Perturbation Method; Analytical Solution.



§ INTRODUCTION

Fragmentation, along with aggregation, represents one of the most basic particulate processes. It refers to the process by which a granule, colloidal cluster, polymer chain or other “particle” breaks into a distribution of fragments. This generic process is encountered in a wide range of natural systems and serves as a model of branching processes in general, in which a population is transformed by replacing a member by a distribution of new members. The physical mechanisms that give rise to breakup or fragmentation (the two terms are used interchangeably in this paper), are generally different between systems but a common mathematical description is possible under the rather general assumption that the rate at which a particle breaks up and the distribution of fragments that are produced depend only on the parent particle and are independent by the presence of other particles in the population. This model views breakup as a first order reaction in the concentration of parent particles, with a distribution of products that depend on the parent particle. Moreover, the physical mechanism that give rise to the aggregation process can be described under the general assumption that the rate at which the particles aggregate and the distribution that are depend on the particles which aggregate together to form a new particle.

The breakage as well as aggregation models will be described by two functions. For defining the breakage model, the rate of breakup of the parent particle, a(m), where m is the size of the particle, and the distribution of fragments, k(m|n) that gives the number of fragments of size m produced by a parent of size n are required. If we take size to refer to the mass of the particle, the fragment distribution obeys the following normalization conditions:


    ∫_0^n k(m|n) dm = f̅,




    ∫_0^n m b(m|n) dm = n.


The first condition gives the average number of fragments, f̅, and is the same for all parent sizes m; the second condition expresses mass conservation between the parent particle and the fragments. A physically realistic fragmentation model requires f̅≥ 2. The governing equation for the size distribution of population that undergoes the above fragmentation process is given by the fragmentation equation,


    ∂ c(m,t)/∂ t =-a(m) c(m,t)_death of particle of size m + ∫_m^∞ a(n) k(m|n) c(n,t) dn_birth of particle of size m,


where c(m,t) is the concentration of particles whose mass is in (m,m+dm) (we normalize the mass m by its mean value at time 0).

Moreover, the aggregation model can be defined with the help of aggregation kernel, g(m,m'), where m and m' are the sizes of the particles which defines the rate at which particles aggregate. The governing equation to track the change in particle size distribution due to aggregation process is defined by aggregation population balance equation,
 
    ∂ c(t,m)/∂ t   =1/2∫_0^m g(m-n, n)c(t,m-n)c(t,n) dn_birth of particle of size m due to aggregation of particles of sizes m-n and n
        -∫_0^∞g(m, n)c(t,m)c(t,n)dn_death of particles of sizes m due to aggregation of particles of sizes m and n,  t∈ [0,T],   x∈(0,∞)

with the initial condition

    c(0,m)=ξ_0(m).

Here c(t,m) represents the concentration of particles of size m at time t. Here, g(m,n) is the aggregation  kernel, which describes the rate at which the particles of sizes m and n coagulate to form a particle of size m+n. It is non-negative and symmetric.

The equations (<ref>) and (<ref>) have been the subject of several investigations and the analytical solutions have been obtained for a small number of special cases. This work, largely theoretical, <cit.> seeks analytical solutions to the fragmentation equation and studies the scaling behavior of the size distribution at long times. A review of these developments and a collection of analytical and solutions is given by Ziff <cit.>.

The analytical solution for aggregation equation can be studied by many researchers for limited number of aggregation kernels. <cit.> derived analytical solution for constant kernel using Taylor polynomials and radial basis functions. Moreover, <cit.> uses the Laplace-Variational iteration to obtain approximate series solution for equation (<ref>). For additive aggregation kernel,  <cit.> proposed analytical solution. However, exact solutions for constant, sum and product kernels are given in <cit.>.







The method of homotopy perturbation has been recently developed to obtain analytical solutions from differential and integral equations and has been successfully used to to solve a a wide range of dynamical problems <cit.>. The purpose of this paper is to formulate the homotopy perturbation method (HPM) in a form appropriate for the fragmentation and aggregation equations and to use it for obtaining analytical solutions. The paper is organized as follows. In section <ref> we introduce the general methodology of HPM. In section <ref>, HPM will be adapt it to the fragmentation and aggregation equations. In next sections <ref> and <ref>, we apply the method to a number of known and new fragmentation and aggregation models. Finally in section <ref> we examine the transient and scaling solution for a  quaternary breakup and aggregation models and relate its behavior.



§ THE HOMOTOPY PERTURBATION METHOD


Recently, the homotopy perturbation method (HPM) has been studied by many researchers for solving linear and nonlinear problems <cit.>. The HPM yields a very rapid convergence of the solution series in most cases, usually only few iterations leading to very accurate solutions.
To illustrate the HPM, we consider the following differential equation

    T(u) - g(r)= 0,     r∈Ω,

with  boundary conditions

    B(u,∂ u/∂ n)= 0,     r∈∂Ω,

where T  is a general differential operator and  B is a boundary operator. Usually the operator T can be decomposed into two parts, a linear operator L and a nonlinear operator N, and expressed as

    L(u)+N(u) - g (r)= 0.

NOTE: If the equation is linear then the non-linear part will be zero, i.e., N(u)=0.

According to HPM, we construct a homotopy  which satisfies

    H[v(r,p)]=(1-p)[L(v(r,p))-L(u_0)]+p[T[v(r,q)] - g(r)]= 0,

where  u_0 is an initial guess to exact solution of (<ref>). When p=0 then L(v(r,0)) = L(u_0)=0,
and when p=1, then T(v(r,1))- g(r) = 0. It is worth noting that as the embedding parameter p increases monotonically from zero to unity, the changing process of p from zero to unity is just that of v(r,p) from u_0(r) to u(r). This is called deformation, and also L(v)-L(u_0) and T(v)- g(r) are called homotopic in topology.

According to the HPM, we can first view the embedding parameter p as a small parameter, and construct the solution as a power series in p, as


    u=∑_k=0^∞ p^kv_k=v_0+p v_1+p^2 v_2+ ⋯.


Substituting (<ref>) in  (<ref>) and then  letting p=1, we obtain the solution as


    f=lim_p→1 u=∑_k=0^∞v_k.


The series  (<ref>) is a convergent for most of the cases  and the rate of convergence depends on the nature of the problem <cit.>.





The idea of the HPM for different types of models <cit.> will be demonstrated in this section. Firstly, the problem linear integral equation will be solved using HPM. Secondly, non-linear integral equation will be taken into account.


Model 1: Consider the linear integral equation:

    c(t)=a(t)+∫_q^t K(u,t,c(u))du.

Here K(u,t,c(u)) should be linear function, i.e., K(u,t,c(u))=∑_j=1^nK_j(u,t)c_i,j(u). For solving above equation using HPM method, a homotopy will be constructed

    (1-p)(C(t)-c_0(t))+p(C(t)-a(t)-∫_q^t K(u,t)C(u)du)=0.

Suppose above equation has the solution of the following form:

    C(t)=C_0(t)+pC_1(t)+p^2 C_2(t)+,

where C_i(t), i=1,2,3, are the functions which will be determined and C_0(t) is the initial approximation to the solution which is taken to be a(t). On substituting (<ref>) into (<ref>) and equating the like powers of p, we have

    p^0    :  C_0(t)-c_0(t)=0,
    
    p^1    :  C_1(t)+c_0(t)-a(t)-∫_q^t ∑_j=1^n K_j(u,t) C_0(u)du=0
    
    p^2    :  C_2(t)-∫_q^t ∑_j=1^n K_j(u,t)C_0(u)du=0
    ⋮


The approximated solution of (<ref>) can be obtained by setting p=1,

    c_i(t)=lim_p → 1 C(t)  = ∑_j=0^∞ C_j(t).


Model 2: Consider the nonlinear  integral equation

    u(x)=g(x)+∫_a^b k(x,s) f(s,u(s)) ds,  x∈Ω.

To apply the HPM, we write (<ref>) as

    L(u)=u(x)-g(x)-∫_a^b k(x,s) f(s,u(s)) ds=0,

with solution y(x) and we construct the homotopy H(u, p) as

    H(u,0) = F(u),    H(u,1)=L(u),

where F(u) is a functional operator with known solution v_0. We choose a convex homotopy by

    H(u,p) = (1-p)F(u)+pL(u)= 0,   p∈ [0,1]

and continuously trace an implicitly defined curve from a starting point H(v_0, 0) to
the solution function H(f, 1). The embedding parameter p increases monotonically from 0 to 1 as a trivial problem F(u)=0 is continuously deformed to the original problem L(u)= 0.

According to the HPM, we can first view the embedding parameter p as a small parameter, and construct the solution as a power series in p, as

    u=∑_k=0^∞ p^kv_k.

 Substituting (<ref>) in  (<ref>) and then  letting p=1, we obtain the solution as

    y=lim_p→1 u=∑_k=0^∞v_k.

The series  (<ref>) is a convergent for most of the cases  and the rate of convergence depends on the nature of the problem <cit.>.

Substitute (<ref>) into (<ref>) and equate the terms with identical powers of p, we obtain

    u_0=g(x),
       u_k+1=∫_a^b k(x,s) H_n ds,   k=0,1,2…




§ IMPLEMENTATION OF HPM


In this section, HPM will be implemented to the aggregation and fragmentation equations. It can be noticed clearly that the breakage equation is linear whereas aggregation equation is non-linear. Therefore, for solving breakage equation, HPM for linear equation will be implemented and for solving aggregation equation, HPM for non-linear equation will be implemented.



 §.§ Implementation of HPM for Breakage Equation

To apply HPM to fragmentation we first express the population balance equation in the integral form,


    ∂ c(t,m)/∂ t =
           - a(m)c(t,m)
           + ∫_m^∞ k(m|n)a(n) c(t,n) dn,


where c_0(m) = c(t,m) is the size distribution at time zero. Next we introduce a new function u=u(t,m; p) of time, size, p and define the homotopy of Eq. (<ref>) as follows:


    (1-p) (∂ c(t,m)/∂ t-∂ c_0(t,m)/∂ t) +
        p{∂ c(t,m)/∂ t
           +a(m)c(t,m) + .
           
    .
           -∫_m^∞ k(m|n)a(n) c(t,n)dn
           } = 0 .


With p=0, Eq. (<ref>) gives u(t,m,0)=c(0,m); with p=1 it gives u(t,m; 1) = c(t,m). Thus by continuously varying p from 0 to 1 we obtain a continuous transformation of u(t,m;p) from the initial distribution c_0 to the actual distribution c(m,t) at time t.

Following HPM, we express u(t,m;p) as a power series in p,


    u(t,m;p) = ∑_k=0^∞ u_k(t,m) p^k,


where the coefficients u_k are function of time and size to be determined. According to the homotopy in Eq. (<ref>) the solution to the fragmentation equation is obtained by setting p=1 in the above series:


    c(m,t) = lim_p→ 1 u(t,m;p) = ∑_k=0^∞ u_k(t,m).


The series is convergent and the rate of convergence depends on the nature of Eq. (<ref>). To obtain the coefficients u_k we insert the series  (<ref>) into (<ref>):


    (1-p)(∂/∂ t∑_k=0^∞ p^k u_k- ∂ c_0/∂ t)
        + p
        {∂/∂ t∑_k=0^∞ p^k u_k+
            a(m)(∑_k=0^∞ p^k u_k)
           . -
           
    .
           ∫_m^∞ k(m|n) a(n)
           (∑_k=0^∞ p^k u_k) dn
        }
        = 0  ,


collect terms in powers of p and set their coefficients to zero. We find:


    k=0           u_0 = c(0,m) ,
    
       k=1         ∂ u_1/∂ t =
          -a(m) u_0(t,m) + ∫_m^∞ k(m|n) a(n) u_0(t,n) dn, 
    ⋮              
    
       k             ∂ u_k/∂ t=
       -a(m) u_k-1(t,m) + ∫_m^∞ k(m|n) a(n) u_k-1(t,n) dn.


The first term, u_0, is equal to the initial distribution and each subsequent term involves the coefficient of previous order. Thus we may obtain these coefficients recursively. The nth order approximation of the solution is obtained by truncating the series past the nth term. Here, however, we are interested in cases for which the infinite series can be expressed in closed analytical form. First we reproduce known solutions to the fragmentation equation using HPM, and we also present new results for a case that has not been solved previously.





 §.§ Implementation of HPM for Aggregation Equation

The homotopy of equation (<ref>) is constructed which satisfies


    H_p   =(1-p)(∂ u(t,m)/∂ t-∂ u_0(t,m)/∂ t)
       +p(∂ u(t,m)/∂ t-1/2∫_0^m a(m-n, n)u(t,m-n)u(t,n) dn+ ∫_0^∞a(m, n)u(t,m)u(t,n)dn )=0



    H_p   =∂ u(t,m)/∂ t-(1-p)∂ u_0(t,m)/∂ t
       +p(-1/2∫_0^m g(m-n, n)u(t,m-n)u(t,n) dn+ ∫_0^∞a(m, n)u(t,m)u(t,n)dn )=0


Inserting u=∑_k=0^∞ p^ku_k in the above equation,

    ∂∑_k=0^∞ p^ku_k(t,m)/∂ t-(1-p)∂ u_0(t,m)/∂ t-p/2∫_0^m g(m-n, n)∑_k=0^∞ p^ku_k(t,m-n)∑_k=0^∞ p^ku_k(t,n) dn
       +p∫_0^∞g(m, n)∑_k=0^∞ p^ku_k(t,m)∑_k=0^∞ p^ku_k(t,n)dn=0


Equate the terms with identical powers of p, we obtain the expression of HPM for aggregation equation:

    ∂ u_1(t,m)/∂ t-∂ u_0(t,m)/∂ t-1/2∫_0^m g(m-n, n) u_0(t,m-n)u_0(t,n) dn
       +∫_0^∞g(m, n)u_0(t,m)u_0(t,n)dn=0
       ∂ u_2(t,m)/∂ t-1/2∫_0^m g(m-n, n) [u_0(t,m-n)u_1(t,n)+u_1(t,m-n)u_0(t,n)] dn
       +∫_0^∞g(m, n)[u_0(t,m)u_1(t,n)+u_1(t,m)u_0(t,n)]dn=0


    ∂ u_3(t,m)/∂ t-1/2∫_0^m g(m-n, n) [u_0(t,m-n)u_2(t,n)+u_2(t,m-n)u_0(t,n)+u_1(t,m-n)u_1(t,n)] dn
       +∫_0^∞g(m, n)[u_2(t,m)u_0(t,n)+u_0(t,m)u_2(t,n)+u_1(t,m)u_1(t,n)]dn=0

⋮
    ∂ u_k(t,m)/∂ t-1/2∫_0^m g(m-n,n) [u_0(t,m-n)u_k-1(t,n)+u_1(t,m-n)u_k-2(t,n)++u_1(t,m-n)
       u_1(t,n)] dn+∫_0^∞g(m, n)[u_k-1(t,m)u_0(t,n)+u_k-2(t,m)u_1(t,n)+
    +u_1(t,m)u_1(t,n)]dn=0

The nth order approximation to the  solution is defined  as 
    ϕ_n=∑_k=0^n u_k.

Note that as n→∞, we have the limit ϕ=∑_k=0^∞ u_k, which is an analytical solution of  (<ref>).
In the next section <ref>, the HPM will used to obtain analytical solution of breakage equation for various breakage kernels with different initial conditions.
































§ APPLICATION OF THE HPM IN FRAGMENTATION EQUATION




 §.§ Case I: Random binary fragmentation with linear selection function


In random binary fragmentation the distribution of fragments is given by


    k(m|n) = 2/n .


Here we obtain the solution with linear fragmentation rate, a(m) = m.  The general form of the factors u_k(m,t) is given by Eq. (<ref>), which now becomes


    ∂ u_k(t,m)/∂ t =
         x u_k-1(t,m) - 2∫_0^∞u_k-1(t,n)/n dn.


Below we obtain the solution for monodisperse and exponential initial distribution.




  §.§.§ Monodisperse initial conditions


We start with u_0=f_0=δ(m-1) and obtain the functions u_k recursively:


    u_0 = δ(m) 
        u_1(m,t) = 2 t θ (1-m)-t x δ (m-1) 
       ⋮
        u_k(m,t) = t^2 (1-m) (-t m)^k-2/(k-2)!θ (1-m),


where θ is the Heaviside step function. The distribution is given by the infinite series


    c(m,t) = ∑_k=0^∞t^2 (1-m) (-t m)^k-2/(k-2)!θ (1-m),


As shown in the Appendix, this is equal to


    c(m,t) =  e^-tm(δ(m-1)+ θ(a-m)(2t+t^2(1-m)) )


This is the same as the result given by Ziff and McGrady <cit.>.




  §.§.§ Exponential initial condition


Starting with u_0(t,m) = f_0(t,m)=e^-m, the u_k are


    u_0(m,t) = e^-m, 
        u_1(m,t) = t (-e^-m) (m-2),
       ⋮
        u_k(m,t) =(-m)^k t^k-2(-2 k t+(k-1) k+t^2) /k! e^-t   .


The size distribution in this case is given by


    c(m,t) = ∑_k=0^∞(-m)^k t^k-2(-2 k t+(k-1) k+t^2) /k! e^-t→ (1+t)^2e^-m(1+t) .


This is the same as the solution given by Ziff and McGrady <cit.>.




 §.§ Case II: Random binary fragmentation with quadratic selection function


In this model the breakup rate is a quadratic function of size, a(m)=m^2 and the fragment size distribution is again given by g(m|n)=2/n and the initial condition is exponential. The general recursion for u_k(m,t) is


    ∂ u_k(t,m)/∂ t = ∫_0^t
         m^2 u_k-1(t,m) - 2∫_0^∞u_k-1(t,n)/n^2 dn.


We obtain solutions for monodisperse and exponential initial conditions:



  §.§.§ Monodisperse initial condition


With u_0(m,t)=δ(m-1) the functions u_k are:


    u_0(m,t) = δ(m-1) 
        u_1(m,t) = 2 t θ (1-m)-t x^2 δ (m-1)
       ⋮
        u_k(m,t) = δ (t-1) (-t^2 x)^k/k!+2 m θ
       (1-t) (-t^2 m)^k-1/(k-1)!


The size distribution is


    c(m,t) = ∑_k=0^∞{δ (t-1) (-t^2 x)^k/k!+2 x θ
                 (1-t) (-t^2 x)^k-1/(k-1)!}
    → e^-tm^2(δ(m-1)+2at θ(1-m)) .


This again recovers the solution of Ziff and McGrady <cit.>




  §.§.§ Exponential initial condition


Starting with c_0(m) = e^-m and solving recursively we obtain


    u_k(m,t) = (-t)^k m^2 k-2k!(-2 k m-2 k+m^2) e^-m .


The solution is constructed by computing the infinite series,


    c(m,t) =∑_k=0^∞(-t)^k m^2 k-2k!(-2 k m-2 k+m^2) e^-m
    → e^-tm^2(δ(m-1)+2 t θ(1-m))  .


This is the same as the result obtained by Ziff and McGrady <cit.>




 §.§ Case III: k-nary fragmentation with power-law rate


In this model the selection function is a(m) = m^α and the fragment size distribution is given by


    k(m|n) = α/n(m/n)^α-2


and number of fragments


    f̅ = α/α -1 .


By varying α between 1 and 2 the number of fragments ranges from ∞ to 2. With α = 2, in particular, this model defaults to the random distribution of fragments with quadratic breakage rate. Equation (<ref>) in this model now is


    ∂ u_k(m,t)/∂ t = x^α u_k-1(t,m)
           -  α∫_m^∞  n m^α-2 u_k-1(t,n) dn.


Using u_0(m,t) = e^-m the functions u_k are


    u_k(m,t) = (-t)^k m^α  k-2k!(-α k-α  k m+m^2) e^-m .


Inserting this result into the infinite series in Eq. (<ref>) we obtain


    c(t,m) = ∑_k=0^∞(-t)^k m^α  k-2k!(-α  k-α  k m+m^2) e^-m
    →
        e^-m-t m^α (1+ α t (m^α-2+m^α-1)) ,


which gives the size distribution at all times. Now, HPM will be apply to solve various models of aggregation equation.





§ APPLICATION OF THE HPM IN AGGREGATION EQUATION





  §.§.§ Constant aggregation kernel g(m, n)=1

The distribution of aggregates for constant aggregation kernel is

    g(m, n)=1.

The general form of the factors u_k(m, t) is given by Eq. (<ref>), which now becomes

    ∂ u_k(t,m)/∂ t-1/2∫_0^m g(m-n, n) [u_0(t,m-n)u_k-1(t,n)+u_1(t,m-n)u_k-2(t,n)++u_1(t,m-n)
        u_1(t,n)] dn+∫_0^∞g(m, n)[u_k-1(t,m)u_0(t,n)+u_k-2(t,m)u_1(t,n)++u_1(t,m)u_1(t,n)]dn=0
    ,
       k=1,2,3...

The solution will be obtained for exponential initial condition, i.e., u_0(m)=e^-m. So, let us initialize the solution with u_0 =e^-m and obtain the functions u_k recursively

    u_0   =e^-m,
    
    u_1   = t/2^1 1!  (m-2)e^-m,
    
    u_2   =t^2/ 2^2 2!(m^2-6 m+6) e^-m ,
    
    ⋮
    
    u_k   = e^-m(k+1)!  t^k/ 2^k∑_r=0^k(-1)^r  m^k-r/r!(k-r+1)!(k-r)!

The partial sum of the series solution is defined   obtained as

    u(t,m)≈ϕ_n(t,m)=∑_k=0^n u_k.

 As n⟶∞, we obtain
 
    u(t,m)=e^-m∑ _k=0^∞t^k /2^k k! HypergeometricU[-k,2,m]

The  exact solution is
u(t,m)=n^2(t)e^-n(t)m, where n(t)=2/t+2 as  in <cit.>.



 §.§ Sum aggregation kernel g(m, n)=m+n

The distribution of aggregates for sum aggregation kernel is defined as

    g(m, n)=m+n

The analytical  solution can be found in <cit.> which is 
    u(t,m)= (1-τ ) e^-(1+τ)m/m√(τ)I[2m√(τ)],
 where
τ =1-e^-t and  I[m] is the modified Bessel function of the first kind.
The general form of the factors u_k(m, t) is given by Eq. (<ref>), which now becomes

    ∂ u_k/∂ t=1/2∫_0^m m    H^1_k-1 dn-∫_0^∞ (m +n)   H^2_k-1dy,    k=1,2,3...




  §.§.§ Exponential initial condition u_0(m)=e^-m

We initialize with u_0 = e^-m and obtain the functions u_k recursively

    u_0   = e^-m,
    
    u_1   =1/2 e^-m t (m^2-2 m-2),
    
    u_2   =1/12 e^-m t^2 (6+18 m-3 m^2-6 m^3+m^4),
    
    u_3   =1/144 e^-m t^3 (-24-168 m-60 m^2+120 m^3+12 m^4-12 m^5+m^6)
    
    u_4   =1/2880 e^-m t^3 (120+1800 m+2100 m^2-1800 m^3-1180 m^4+360 m^5+70 m^6-20 m^7+ m^8)
    ⋮

It is difficult to find the generalized term for the solution of sum kernel. So, we approximate the solution by taking the summation of some terms. The comparison of exact solution with the approximate solution is represented in Figure <ref>. The comparison reveals that by considering the sum of four terms in the series, exact solution can be approached.







 §.§ Product aggregation kernel g(m, n)=mn

The distribution of aggregates for product aggregation kernel is

    g(m, n)=mn

The general form of the factors u_k(m, t) is given by Eq. (<ref>), which now becomes

    ∂ u_k(t,m)/∂ t-1/2∫_0^m (m-n)n [u_0(t,m-n)u_k-1(t,n)+u_1(t,m-n)u_k-2(t,n)++u_1(t,m-n)
        u_1(t,n)] dn+∫_0^∞mn[u_k-1(t,m)u_0(t,n)+u_k-2(t,m)u_1(t,n)++u_1(t,m)u_1(t,n)]dn=0
    ,
       k=1,2,3...

The solution will be obtained for two different initial conditions, i.e., for u_0(m)= e^-m and u_0(m)= e^-m/m.



  §.§.§ Exponential initial condition, u_0(m)= e^-m

Its analytic solution follows <cit.> as  
    u(t,m)=∑ _k=0^∞ e^-(t+1)mt^k m^3k/(k+1)!  Γ(2k+2)








We start with u_0 =e^-m and obtain the functions u_k recursively

    u_0   = e^-m,
    
    u_1   =(tm)^1/2! 3!  (-12+m^2)e^-m,
    
    u_2   =(tm)^2/ 3!5! (360-60 m^2+m^4)e^-m,
    
    ⋮
    
    u_k   =  2e^-m(mt)^k ∑_r=0^k(-1)^r  m^2k-2r/r!(k-r)!  (2k-2r+2)!

The partial sum of the series solution is defined   obtained as

    u(t,m)≈ϕ_n(t,m)=∑_k=0^n u_k.

Now, the solution for another initial condition.


  §.§.§ When u_0(m)= e^-m/m

The analytical solution can be found in <cit.>

    u(t,m)=I[2m√(t)]/m^2√(t)  e^-Tm,

where
T = {[      1+t ,       t≤1,;     2√(t), otherwise. ].


Here I is the modified Bessel's first kind function:

    I[m]=1/π∫_0^πe^(m cosθ)cosθ dθ.

We initialize with u_0 = e^-m/m and obtain the functions u_k recursively

    u_0   =e^-mm^-1/1,
    
    u_1   =t/2 (1!)^2  (m-2)e^-m ,
    
    u_2   =t^2 m/3 (2!)^2(m^2-6 m+6)e^-m ,
    
    ⋮
    
    u_k   = e^-m x^k-1 t^k ∑_r=0^k(-1)^r m^k-r/r!(k-r+1)!(k-r)!.

The partial sum of the series solution is obtained as

    u(t,m)≈ϕ_n(t,m)=∑_k=0^n u_k.

 As n⟶∞, we obtain

    e^-m∑ _k=0^∞t^k m^k-1/k! (k+1)! HypergeometricU[-k,2,m]





§ DISCUSSION


Cases 1 and 2, which describe random binary breakup with linear and quadratic rate functions essentially form the “classical” model in fragmentation. Their solutions have been obtained previously in the literature multiple times by various methods, including probabilistic treatment <cit.>, and solution of the equivalent discrete problem <cit.>. The kernel of Case 3 is a special form of a kernel whose general form is


    k(m|n) =
        (1-ψ ) λm^λ -1/n^λ +
       ψγm^γ -1/n^γ .


This kernel was introduced by Austin <cit.> to model grinding. In cumulative form it is the sum of two power-law distributions, of which one is the dominant distribution for large m (Austin gave this kernel with different coefficients but the two forms can be shown to be the same under appropriate relationships among the coefficients). The same kernel has been the subject of theoretical investigations of the scaling behavior of the fragmentation equation <cit.> and of the so-called shattering transition <cit.>, a behavior analogous but opposite to gelation that is manifested by the presence of a finite fraction of the mass of the system in particles of zero size.  Ziff and McGrady <cit.> gave the solution of the Austin kernel with power-law breakage rate for monodisperse initial conditions. Oukouomi Noutchie and Doungmo Goufo obtained the same solution by Laplace transform <cit.>.



In principle, given the linearity of the fragmentation equation the solution for any other initial condition can be reconstructed from the solution to the monodisperse problem. The monodisperse solution, however, employs the confluent hypergeometric function and its manipulation is not trivial. The solution obtained here in Eq. (<ref>) for exponential conditions is compact and represents a new result. This solution obeys the scaling derived by Ziff and McGrady <cit.>, as can be easily demonstrated. Dividing Eq. (<ref>) by m^2/α and setting m=(z/t)^1/α we find


    α c(m,t)/m^2 =
        1/α{
        t^-2/α e^t^-1/α(
           - z^1/α)-z(α  t^1/α z^(α+1)/α
           + α  z t^2/α
           + z^2/α)
           } ,


whose limit at long t is


    α c(m,t)/m^2→ z e^-z≡Φ(z) ,


and Φ(z) is the scaling function. The scaling limit of the distribution is


    c(m,t) ∼α m^2Φ(t m^a) .







The scaling function given by Ziff and McGrady is Φ(z) = z^γ/α e^-zΓ(γ/α) (the Filippov model in Table 1 of Ref. <cit.> ), which for γ=α reduces to the result obtained here.

To illustrate these results we examine the solution for α = 4/3. The fragmentation rate is a(m)=m^4/3 and the distribution of fragments is


    k(m|n) = 4/3 m^2/3 n^3/2,


with an average number of fragments f=4. This distribution is very asymmetric with a large fraction of the population in the small size region (Fig. <ref>). There is no shattering, however, because the fragmentation rate decreases as the size approaches zero (indeed our solution does include shattering as the condition f≥ 2 constrains α to be positive). The size distribution is shown in Fig. <ref>. After an initial transition the distribution in log-log coordinates reaches an invariant form that is simply translated to lower average sizes. This is a manifestation of the scaling limit.

Finally, to examine the convergence of the series we calculate the solution by retaining only a finite number of terms k. Figure <ref> compares truncated approximations of the full solutions for α=3/2 at t=10 and k_max=2,3,12 and 13 (k_max is the maximum order of the term retained in the series). The general feature of the truncation is that the approximation is excellent for 0<m<m_max, where m_max is the size beyond which the approximation breaks down and is pushed to larger sizes as the order of the truncation increases. Even-order truncations provide stable approximations approximations of the size distribution and it is interesting to point out that in this particular case, even a very low order approximation such as k_max=2 represents the distribution very well over two decades in c(m).






















































































§ CONCLUSIONS


The homotopy perturbation method has been successfully applied to the fragmentation as well as aggregation population balance equation. The HPM provides a systematic methodology to solve for the size distribution and was shown to reproduce all known solutions and to lead to new results for a special case of the Austin kernel and also reproduce analytical solution for various models of fragmentation and aggregation.





apsrev
