**2. Materials and Methods**

Suppose *ud* represents the number density function, then a general PBE is given by [1,2]:

$$\frac{\partial u\_d(T, V\_p)}{\partial T} = -\underbrace{\frac{\partial [G(T, V\_p)u\_d(T, V\_p)]}{\partial V\_p}}\_{\text{Convoth-term}} + \underbrace{\underbrace{Q\_{\text{muc}}(T, V\_p)}\_{\text{Nuckation-term}} + \underbrace{Q\_{\text{AgS}}^\pm(T, V\_p)}\_{\text{Agygness-term}}}\_{\text{Agygness-term}} + \underbrace{Q\_{\text{Broak}}^\pm(T, V\_p)}\_{\text{Broakage-term}} + \underbrace{Q\_{\text{diss}}(T, V\_p)}\_{\text{Dissolution-term}} \tag{1}$$

where R+ = (0, <sup>∞</sup>). The variable *T* represents the time and *Vp* may be size, length, or composition. In this article, it represents the particle volume. In the above equation, each term has its specific definition. These terms are given by

$$Q\_{\rm agg}^{\pm}(T, V\_p) = B\_{\rm agg}(T, V\_p) - D\_{\rm agg}(T, V\_p) \tag{2}$$

where

$$B\_{\rm agg}(T, V\_p) = \frac{1}{2} \int\_0^{V\_p} \beta \Big(T\_\prime V\_p - V\_{p\prime}' V\_p'\Big) \mu\_d \Big(T\_\prime V\_p - V\_p'\Big) \mu\_d \Big(T\_\prime V\_p'\Big) dV\_p'$$

and

$$D\_{\mathrm{agg}}(T, V\_p) = \bigcup\_{0=1}^{\infty} \beta \big(T, V\_{p\prime}, V\_p' \big) \mu\_d \big(T, V\_p \big) \mu\_d \big(T, V\_p' \big) dV\_p'$$

Here, *<sup>B</sup>*agg*<sup>T</sup>*, *Vp* represents the birth of particles of volume *Vp* resulting for amalgamation of two particles with respective volume *Vp* and *Vp* − *<sup>V</sup>p*, where *Vp* ∈ 0, *Vp* and *<sup>D</sup>*agg*<sup>T</sup>*, *Vp* is the death of particles, describing the decrease in particle volume *Vp* by aggregation with other particles of any volume. The aggregation kernel <sup>β</sup>*<sup>T</sup>*, *Vp*, *Vp* is the rate at which aggregation of two particles *Vp* and *Vp* produces a particle volume *Vp* + *<sup>V</sup>p*.

$$Q\_{\rm break}^{\pm}(T, V\_p) = B\_{\rm break}(T, V\_p) - D\_{\rm break}(T, V\_p) \tag{3}$$

where

$$B\_{\rm break}(T, V\_p) = \int\_{V\_p}^{\infty} b\left(T, V\_{p^\*}, V\_p'\right) \mathcal{S}\left(V\_p'\right) \mu\_d\left(T, V\_p'\right) dV\_p'$$

and

$$D\_{\text{break}}(T, V\_p) = S(V\_p) \mu\_d(T, V\_p)$$

Here, *<sup>B</sup>*break*<sup>T</sup>*, *Vp* represents the birth of new particles during the breakage process and the breakage function *<sup>b</sup><sup>T</sup>*, *Vp*, *Vp* is the probability density function for the formation of particle volume *Vp* from particle volume *<sup>V</sup>p*, whereas, *<sup>D</sup>*break*<sup>T</sup>*, *Vp* is the death of particles, and *<sup>S</sup><sup>V</sup>p* is the selection function describing the rate at which the particles are selected to break.

$$Q\_{\rm diss}(T, V\_p) = \frac{\dot{V}}{V\_c} h(V\_p) \mu\_d(T, V\_p) \tag{4}$$

*<sup>Q</sup>*diss*<sup>T</sup>*, *Vp* represents the dissolution term, *Vc* is the volume of crystallizer, .*V* is the volumetric flow rate from the crystallizer to dissolution unit, *hVp*is the dissolution function describing dissolution of small particles below some critical size. The population balance equation can be simplified through introducing the moment function. The *kth* moment of the population density function is mathematically written in the form:

$$m\_k(T) = \bigcup\_{0}^{\infty} V\_p^k \mu\_d(T, V\_p) dV\_p \tag{5}$$

The first and second moments *<sup>m</sup>*0(*T*), *<sup>m</sup>*1(*T*) denote particle population and volume, respectively, at any instant *T*. Multiply *Vkp* to left- and right-hand side of Equation (1) and at that moment integrate it over the volume *Vp*, so we obtain the following equation:

$$\frac{d\eta\_1(T)}{dT} = \overline{Q\_{\text{growth}}}(T, V\_p) + \overline{Q\_{\text{nucleation}}}(T, V\_p) + \overline{Q\_{\text{approx}}}(T, V\_p) + \overline{Q\_{\text{twoakage}}}(T, V\_p) + \overline{Q\_{\text{twoakage}}}(T, V\_p) + \overline{Q\_{\text{dissulation}}}(T, V\_p) \tag{6}$$

where

$$\overline{G\_{\text{growth}}}(T, V\_p) = \int\_0^\infty V\_p^k \frac{\partial [G(T, V\_p) u\_d(T, V\_p)]}{\partial V\_p} dV\_p$$

$$\overline{Q\_{\text{nucleation}}}(T, V\_p) = \int\_0^\infty V\_p^k Q\_{\text{nuc}}(T, V\_p) dV\_p$$

$$\overline{Q\_{\text{aggregation}}}(T, V\_p) = \overline{B\_{\text{agg}}}(T, V\_p) - \overline{D\_{\text{agg}}}(T, V\_p)$$

The aggregation terms *<sup>B</sup>*agg(*<sup>T</sup>*, *Vp*) and *<sup>D</sup>*agg(*<sup>T</sup>*, *Vp*) are mathematically defined as

$$
\overline{B\_{\rm agg}}(T, V\_p) = \frac{1}{2} \int\_0^\infty V\_p^k \int\_0^{V\_p} \beta \left(T, V\_p - V\_{p'}', V\_p'\right) \mu\_d \left(T, V\_p - V\_p'\right) \mu\_d \left(T, V\_p'\right) dV\_p dV\_p
$$

$$
\overline{D\_{\rm agg}}(T, V\_p) = \int\_0^\infty V\_p^k \mu\_d \left(T, V\_p\right) \int\_0^\infty \beta \left(T, V\_{p'}, V\_p'\right) \mu\_d \left(T, V\_p'\right) dV\_p' dV\_p
$$

$$
\overline{Q\_{\rm breakage}}(T, V\_p) = \overline{B\_{\rm break}}(T, V\_p) - \overline{D\_{\rm break}}(T, V\_p)
$$

*Processes* **2019**, *7*, 453

whereas *<sup>B</sup>*break(*<sup>T</sup>*, *Vp*) and *<sup>D</sup>*break(*<sup>T</sup>*, *Vp*) are birth and death functions because of breakage term and are given by

$$
\overline{B\_{\rm break}}(T, V\_p) = \bigcap\_{0}^{\infty} V\_p^k \bigcap\_{V\_p}^{\infty} b\{T, V\_p, V\_p'\} S(V\_p') \mu\_d(T, V\_p') dV\_p V\_p
$$

$$
\overline{B\_{\rm break}}(T, V\_p) = \int\_0^{\infty} V\_p^k S(V\_p) \mu\_d(T, V\_p) dV\_p
$$

$$
\overline{Q\_{\rm dissolution}}(T, V\_p) = \frac{\dot{V}}{V\_c} \int\_0^{\infty} V\_p^k h(V\_p) \mu\_d(T, V\_p) dV\_p
$$

Due to aggregation, the upper limit in the birth function is not infinity; therefore, we cannot apply quadrature method of moment for this integral. To apply the QMOM, we have to convert the upper limit to infinity. Here, we introduce a Heaviside step function *H Vp* − *V p* to solve the integral such that *H Vp* − *V p* = 0 when *Vp* − *V p* < 0 and *H Vp* − *V p* = 1 otherwise. As a result, we will find the limits of integration over *V p* from 0, *Vp* to (0, <sup>∞</sup>). Thus, by applying it to the birth function, we obtain the following equation:

$$\begin{split} \frac{1}{2} \int\_{0}^{\infty} V\_{p}^{k} dV\_{p} \int\_{0}^{\infty} \mathbf{H} \Big( V\_{p} - V\_{p}^{\prime} \Big) \boldsymbol{\beta} \Big( T, V\_{p} - V\_{p}^{\prime}, V\_{p}^{\prime} \Big) \boldsymbol{\mu}\_{d} \Big( T, V\_{p} - V\_{p}^{\prime} \Big) \boldsymbol{\mu}\_{d} \Big( T, V\_{p}^{\prime} \Big) dV\_{p} \\ = \frac{1}{2} \int\_{0}^{\infty} \boldsymbol{\mu}\_{d} \Big( T, V\_{p}^{\prime} \Big) dV\_{p}^{\prime} \int\_{0}^{\infty} \Big( \boldsymbol{\mu} + V\_{p}^{\prime} \Big)^{k} \Big{\beta} \Big( T, \boldsymbol{\mu}, V\_{p}^{\prime} \Big) \boldsymbol{\mu}\_{d} \Big( T, \boldsymbol{\mu} \Big) dV\_{p} \end{split} \tag{7}$$

On the right side of Equation (7), we have switched the order of integration and made the replacement *u* = *Vp* − *V p*. Lastly, we have replaced *Vp* for *u* with no damage of generality and caught the necessary outcome for birth due to aggregation which is given by

$$\overline{B\_{\rm agg}}(T, V\_p) = \frac{1}{2} \int\_0^\infty u\_d(T, V'\_p) \int\_0^\infty (V\_p + V'\_p) \left\{(T, V\_{p'}, V'\_p) u\_d(T, V\_{p'}) dV\_p dV'\_{p'} \right\}$$

After substituting all of the above terms in Equation (6), we ge<sup>t</sup> the system of di fferential equations:

$$\begin{array}{cl} \frac{d\boldsymbol{u}\_{d}(T)}{dT} &=& \int\_{0}^{\infty} kV\_{p}^{k-1} \mathbb{C}\{T, V\_{p}\} \boldsymbol{u}\_{d}\{T, V\_{p}\} dV\_{p} + \int\_{0}^{\infty} V\_{p}^{k} \mathbb{Q}\_{\text{nuc}}\{T, V\_{p}\} dV\_{p} \\ &+ \frac{1}{2} \int\_{0}^{\infty} \boldsymbol{u}\_{d}\{T, V\_{p}^{\prime}\} \begin{pmatrix} \boldsymbol{V}\_{p} + \boldsymbol{V}\_{p} \end{pmatrix} \boldsymbol{\beta}\{T, V\_{p}, V\_{p}^{\prime}\} \boldsymbol{u}\_{d}\{T, V\_{p}\} dV\_{p} dV\_{p} \\ &\overset{\text{or}}{\quad} - \int\_{0}^{\infty} \boldsymbol{V}\_{p}^{k} \boldsymbol{u}\_{d}\{T, V\_{p}\} \int\_{0}^{\infty} \boldsymbol{\beta}\{T, V\_{p}, V\_{p}^{\prime}\} \boldsymbol{u}\_{d}\{T, V\_{p}^{\prime}\} dV\_{p} dV\_{p} \\ &+ \int\_{0}^{\infty} \boldsymbol{V}\_{p}^{k} \int\_{0}^{\infty} \boldsymbol{b}\{T, V\_{p}, V\_{p}\} \boldsymbol{S}\{V\_{p}\} \boldsymbol{u}\_{d}\{T, V\_{p}\} dV\_{p} dV\_{p} \\ &- \int\_{0}^{\infty} \boldsymbol{V}\_{p}^{k} \boldsymbol{S}\{V\_{p}\} \boldsymbol{u}\_{d}\{T, V\_{p}\} dV\_{p} + \underset{V\_{c}}{V\_{c}} \int\_{0}^{\infty} \boldsymbol{V}\_{p}^{k} \boldsymbol{h}\{V\_{p}\} \boldsymbol{u}\_{d}\{T, V\_{p}\} dV\_{p} \end{array} \tag{8}$$

In addition, with solid phase, the liquid phase yields ordinary di fferential equation for the solute mass in the form:

$$\frac{dm(T)}{dT} = \dot{m\_{\rm in}}(T) - \dot{m\_{\rm out}}(T) - 3\sigma\_c k\_{\rm v} \int\_0^\infty V\_p^2 \mathbb{G}(T, V\_p) \mu\_d(T, V\_p) dV\_p \tag{9}$$

with *m*(0) = *m*0 where σ*c* is the density of crystals and *kv* is a volume shape factor. . *<sup>m</sup>*in(*T*) − . *<sup>m</sup>*in(*T*)− is the incoming mass flux from dissolution unit to crystallizer and outgoing mass flux from crystallizer to dissolution unit, respectively, which is mathematically defined by

$$m\_{\rm out}^{\dot{}}(T) = \frac{m(T)}{m(T) + m\_{\rm solv}(T)} \sigma\_{\rm sol}(T) \dot{V} \tag{10}$$

where *<sup>m</sup>*solv(*T*) is the mass of solvent and <sup>σ</sup>*sol*(*T*) is the density of the solution. *Tp* = *Vp V* ≥ 0 is the residence time in the dissolution unit, where *Vp* represents volume of the pipe.

Here, we will consider the Gaussian quadrature method to solve the complicated integral terms appearing in Equations (8) and (9). Consequently, a closed-form system of moments is found. Further, we have accurately and e fficiently solved the closed-form system with an ODE solver. The approximation of definite integral with the help of quadrature rule is an important numerical aspect. For this purpose, first we calculate the function at definite points and then use the formula of weight function which gives approximation of the definite integral. To approximate a definite integral given in Equations (8) and (9), first we find the values of the function at a set of equidistant points. After evaluating, the weight function is multiplied to approximate the integrals. In this rule, there is no limitation of choosing abscissas and weights. This rule also works for the points which are not likewise spread out. Let us assume the integral of the formula *b a* ψ *Vp ud Vp dVp*. We can find the weights *wj* and abscissas *Vkpi*by approximating the definite integral as

$$\int\_{a}^{b} \psi(V\_p)\mu\_d(V\_p)dV\_p = \sum\_{j=1}^{N} w\_j \mu\_d(V\_{p\_j}) \tag{11}$$

where provided it is exact and *ud Vp* is smooth function. A set of orthogonal polynomials is necessary for making certain *nth* order orthogonal polynomial and it should contain only one polynomial of order *j* for *j* = 1, 2, 3, .... We can define the scalar product of the two functions *k Vp* and *h Vp* over a weight function ψ *Vp* by

$$ = \int\_{a}^{b} \psi(V\_P) h(V\_P) k(V\_P) dV\_P \tag{12}$$

If we find the zero scalar product of any two functions *h* and *k*, then these functions are termed as orthogonal functions. It is also noted that supplementary information will be required for obtaining abscissas and weights if the classical weight function ψ *Vp* is not provided. For this purpose, we have used the moment:

$$m\_i(T) = \bigcap\_{0}^{\infty} V\_p^i u\_d(T, V\_p)dV\_p \approx \sum\_{j=1}^{N} V\_{p\_j}^j w\_j \tag{13}$$

In the above equation, *ud T*, *Vp* is used as a weight function ψ *Vp* . *Vi p* given in Equation (13) denotes the polynomial *ud Vp* of *ith* order. To find *n* abscissas weights, we need 2*n* moments from *m*0 to *m*2*n*−1. For *i* ≤ 2*n* − 1 this approximation will be exact. For simplicity and accurate approximations, we set *n* = 3. As a result, six-moment system for *m*0, ...... , *m*5 is calculated. To discuss the procedure, we consider ODEs (8) and (9). Using Equation (13) in Equations (8) and (9), we obtain the following equations:

$$\begin{array}{ll} \frac{dm\_{\mathcal{V}}(T)}{dT} &= k \sum\_{i=1}^{N} \left( V\_{p\_{i}} \right)^{k-1} w\_{i} \mathbf{G} \left( T, \ V\_{p\_{i}} \right) + a\_{\text{max}}^{(k)} \\ &+ \frac{1}{2} \sum\_{i=1}^{N} w\_{i} \sum\_{j=1}^{N} \left( V\_{p\_{i}} + V\_{p\_{j}} \right)^{k} \boldsymbol{\beta} \{ T, V\_{p\_{i}}, V\_{p\_{j}} \} w\_{j} \\ &- \sum\_{i=1}^{N} \left( V\_{p\_{i}} \right)^{k} w\_{i} \sum\_{j=1}^{N} \boldsymbol{\beta} \{ T, V\_{p\_{i}}, V\_{p\_{j}} \} w\_{j} \\ &+ \sum\_{i=1}^{N} \left( V\_{p\_{i}} \right)^{k} w\_{i} \sum\_{j=1}^{N} \boldsymbol{\beta} \{ T, V\_{p\_{i}}, V\_{p\_{j}} \} \mathbf{S} \{ V\_{p\_{j}} \} w\_{j} \\ &+ \frac{\dot{V}}{\nabla\_{c}} \sum\_{i=1}^{N} w\_{i} \left( V\_{p\_{i}} \right)^{k} h \{ V\_{p\_{i}} \} \qquad k = 0, 1, 2, \dots, \dots \end{array} \tag{14}$$

where *a*(*k*) *nuc* = ∞ 0 *VpkQnucdVp* and

$$\frac{dm(T)}{dT} = \dot{m\_{\rm in}}(T) - \dot{m\_{\rm out}}(T) - 3\sigma\_c k\_v \sum\_{i}^{N} w\_i V\_p^2 \mathcal{G}\{T\_\prime, V\_{p\rm i}\} \tag{15}$$

Our next step is computing the quadrature points *Vpi*, *Vpj* and the quadrature weights *wi* and *wj*. The quadrature points *Vpi*, *Vpj* are obtained from the roots of orthogonal polynomials. The construction of orthogonal polynomials is described as follows:

$$p\_{-1} = 0, p\_0 = 1, p\_j = (V\_p - \alpha\_j)p\_{j-1} - \beta\_j p\_{j-1} \tag{16}$$

with

$$\alpha\_{j} = \frac{}{}, j = 1,2,... \tag{17}$$

$$\beta\_{\vec{\beta}} = \frac{}{}, j = 2,3,... \tag{18}$$

Since the function *ud<sup>T</sup>*, *Vp* is used as a weight function <sup>ψ</sup>*Vp*, so from Equation (12) we have

$$ = \int\_a^b u\_d(T\_\prime V\_p) p\_j^2 dV\_p$$

Using all of the above definitions, the orthogonal polynomials of any order can be obtained. Our first step is to find the roots of the *nth* order polynomial. The quadrature points *Vp* are then obtained from these roots. To explain the procedure of deriving these orthogonal polynomials, we have calculated *<sup>p</sup>*<sup>1</sup>*Vp*, *<sup>p</sup>*<sup>2</sup>*Vp*, *<sup>p</sup>*<sup>3</sup>*Vp*. *<sup>p</sup>*1*Vp* as

$$p\_1(V\_p) = (V\_p - \alpha\_1)p\_0 = (V\_p - \alpha\_1)^\*$$

First, α1 will be calculated, which is given by

$$a\_1 = \frac{}{} = \frac{\int\_0^\infty V\_p u\_d(T, V\_p) p\_0^2 dV\_p}{\int\_0^\infty u\_d(T, V\_p) p\_0^2 dV\_p} = \frac{m\_1(T)}{m\_0(T)}$$

$$p\_1(V\_p) = V\_p - \frac{m\_1}{m\_0} \tag{19}$$

so Now

$$p\_2(V\_p) = (V\_p - \alpha\_2)p\_1 - \beta p\_1 \tag{20}$$

Again from Equations (17) and (18) we have

$$\alpha\_{2} = \frac{\,\_{\text{c}}V\_{p}p\_{1}/p\_{1}\,\text{s}}{\,\_{\text{c}}^{\text{in}}/p\_{1}\,\text{s}} = \frac{\int\_{0}^{\infty}V\_{p}u\_{d}(T,V\_{p})p\_{1}^{2}dV\_{p}}{\int\_{0}^{\infty}u\_{d}(T,V\_{p})p\_{1}^{2}dV\_{p}} = \frac{\int\_{0}^{\infty}V\_{p}u\_{d}(T,V\_{p})\left(V\_{p} - \frac{m\_{1}}{m\_{0}}\right)^{2}dV\_{p}}{\int\_{0}^{\infty}u\_{d}(T,V\_{p})\left(V\_{p} - \frac{m\_{1}}{m\_{0}}\right)^{2}dV\_{p}} = \frac{m\_{3}m\_{0}^{2} - 2m\_{0}m\_{1}m\_{2} + m\_{1}^{3}}{m\_{2}m\_{0}^{2} - m\_{0}m\_{1}^{2}}\tag{21}$$

and

$$\beta\_2 = \frac{}{} = \frac{\int\_0^\infty V\_p \mu\_d(T\_\star V\_p) \left(V\_p - \frac{m\_1}{m\_0}\right)^2 dV\_p}{\int\_0^\infty \mu\_d(T\_\star V\_p) dV\_p} = \frac{m\_2 m\_0 - m\_1^2}{m\_0^2} \tag{22}$$

so from Equation (20) we have

$$p\_2(V\_p) = \frac{V\_p^2(m\_0 m\_2 - m\_1^2) + V\_p(m\_1 m\_2 - m\_0 m\_3) + m\_1 m\_3 - m\_2^2}{m\_0 m\_2 - m\_1^2} \tag{23}$$

Proceeding in a similar way, we can calculate the polynomial of higher order. The third-order polynomial *<sup>p</sup>*3*Vp* is given by the following equation:

$$\begin{array}{ll} p\_3(V\_p) &= V\_p^{-3} + \frac{\left(m\_2 m\_4 m\_1 - m\_0 m\_4 m\_3 + m\_2 m\_3 m\_5 + m\_3^2 m\_1 - m\_2^2 m\_5\right) V\_p^2}{m\_2^3 - m\_2 m\_4 m\_0 - 2m\_2 m\_3 m\_1 + m\_3^2 m\_0 + m\_4 m\_1^2} \\ &+ \frac{\left(m\_2 m\_5 m\_1 + m\_0 m\_4^2 - m\_3 m\_5 m\_3 - m\_4 m\_3 m\_1 - m\_2^2 m\_4 + m\_3^2 m\_2\right) V\_p}{m\_2^3 - m\_2 m\_4 m\_0 - 2m\_2 m\_3 m\_1 + m\_3^2 m\_0 + m\_4 m\_1^2} \\ \end{array} \tag{24}$$

The roots of the selected polynomial will give us the abscissas *Vp*. Next, the weights *wi* will be calculated. According to Press et al. [15], the expression for the weight function is given by

$$w\_i = \frac{\langle p\_{N-1}/p\_{N-1} \rangle}{p\_{N-1}\left(x\_j\right)p\_N'\left(x\_j\right)} \quad i = 1, 2, \dots, \dots, N \tag{25}$$

where *N* is the order of the selected polynomial. At last, the resulting system of ODEs is then solved by any standard ODE solver in MatLab.
