# Population Balance Solution Using Fixed Point Convergence Acceleration for Steady-State Reactor Networks

Cristian Camilo Ruiz Vasquez, Noureddine Lebaz, Isabelle Ramière, Murielle Bertrand, Denis Mangin

► To cite this version:

Cristian Camilo Ruiz Vasquez, Noureddine Lebaz, Isabelle Ramière, Murielle Bertrand, Denis Mangin. Population Balance Solution Using Fixed Point Convergence Acceleration for Steady-State Reactor Networks. EPJ Web of Conferences, 2025, 6th International ATALANTE Conference on Nuclear Chemistry for Sustainable Fuel Cycles (ATALANTE-2024), 317, pp.01010. 10.1051/epj-conf/202531701010. cea-04925971

HAL Id: cea-04925971

https://cea.hal.science/cea-04925971v1

Submitted on 3 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# POPULATION BALANCE SOLUTION USING FIXED POINT CONVERGENCE ACCELERATION FOR STEADY-STATE REACTOR NETWORKS

Cristian Camilo Ruiz Vasquez<fnref n="1" />, Noureddine Lebaz<fnref n="2" />, Isabelle Ramière<fnref n="3" />, Murielle Bertrand<fnref n="1" />*, Denis Mangin<fnref n="2" />

**Abstract.** The simulation of industrial crystallizers needs the implementation of multi-compartment models to handle the spatial heterogeneities at lower cost. In this context the study of recycle streams appears as an intermediate step to solve complex reactor networks. In this work, an accelerated fixed-point algorithm is applied to the solution of the crystal size distribution and the liquid phase composition in a structure including two compartments and one recycle stream. The numerical strategy, including fixed-point formulation, acceleration method, initial condition and convergence criteria, is introduced. One of the strengths of the proposed algorithm is that it does not depend on the numerical method used to solve the population balance equation in each compartment. As a case study, the numerical method is performed on the uranium oxalic precipitation. Thanks to this approach, local values of crystallisation kinetics, mass transfer and crystal size distribution are obtained. The evolution of the actinide concentration and the average crystal size follow the theory of the crystallisation mechanisms. Finally, the very satisfactory numerical performances open the way to simulate more complex and demanding reactor network structures.

## 1. Context and numerical issue

In process engineering, most unit operation equipment is characterised by complex hydrodynamics leading to heterogeneous properties all over the geometrical volume. To model real operations, the multi-compartment model presents an alternative between the high resource-consuming CFD simulations and the global homogeneous models [1]. The development of a multi-compartment model can be performed following three successive steps:

* First, the size and number of compartments and the flowrates between the compartments are fixed.

* Second, the physicochemical phenomena taking place in every single homogeneous volume are defined.

* Finally, a suitable solving methodology is performed to consider the interactions between the different homogeneous compartments.

Following this methodology, this research work appears as an intermediate step in the simulation of the oxalic precipitation in a vortex reactor. First, the hydrodynamics of the reactor was studied numerically via Large-Eddy Simulations (LES) [2]. As a result, the complex hydrodynamics taking place in the reactor can be described in a simplified way by five homogeneous compartments (Figure 1c)). In this concern, a robust numerical methodology was developed to solve the population balance equation in each homogenous compartment (Figure 1a)) [3,4].

The current work focuses on developing a numerical methodology to solve the reactor network. To this end, we propose to consider the recycling stream as an elementary block in the solution process and to reuse the algorithm previously introduced to solve each compartment [3,4]. A global fixed-point strategy manages the convergence of the global network solution. In the next section, the fixed-point method used to solve the entire crystal size distribution is described. Then in Section 3, the case study

<footnote n="1">¹CEA, DES, ISEC, DMRC, Université Montpellier, Marcoule, France</footnote>

<footnote n="2">²Université Claude Bernard Lyon 1, CNRS, LAGEPP UMR 5007, Villeurbanne, France</footnote>

<footnote n="3">³CEA, DES, IRESNE, DEC, F- 13108 Saint-Paul Lez Durance, France</footnote>

<footnote>* Corresponding author</footnote>
under consideration is introduced: the uranium oxalate precipitation including nucleation and growth in a simple recycling structure (cf. Figure 1b)).

<figure><img src="image_2.png" /><figcaption>Figure 1. Different reactor networks strategies ($F_i$ and $R_i$ state for the $i^{th}$ stream and the $i^{th}$ reactor respectively).</figcaption></figure>

## 2. Numerical method

The numerical strategy presented in this section aims to determine the limiting reactant concentration $A$ (mol m⁻³) and the crystal size distribution $N$ (m⁻³) of the recycle stream. The internal algorithms used to solve the population balance equation in each compartment ($f_{R_i}$) determine the objective function:

$$ (A, N) = f_{R_n} \left( f_{R_{n-1}} \dots f_{R_1} (A_1, N_1) \right) \qquad (1) $$

Whit $i$ the number of compartments. The reader must note that the dependence on the other process variables such as volume, temperature and flowrate is omitted for clearness purposes.

### 2.1 Fixed point algorithm and acceleration method

When one of the feeding streams is obtained by recycling a part of the product stream (see Figure 1b) for example), at the steady-state the input variables of the entire system are partially set by the product characteristics itself. From a mathematical point of view, the self-dependence of an output variable can be expressed as a fixed-point problem: $X = f(X)$, with $X$ the unknown vector. In the precipitation context the fixed point problem becomes: $(A, N) = (A, N)$. In the particular case of Figure 1b), we hence have:

$$ (A_5, N_5) = f(A_5, N_5) = f_{R_2} \left( f_{R_1} (A_5, N_5, A_1, N_1) \right) \qquad (2) $$

In this case, the outgoing crystal size distribution $N_5$ and the concentration of the limiting reagent ($A_5$, namely the iterating variable) are defined by the solution of the population balance equations of reactor $R_1$ followed by reactor $R_2$ ($f_{R_2}(f_{R_1})$). These population balance equations depend on the crystal size distribution $N_5$ and the concentration in the feed stream $A_5$ itself.

A well-known numerical method to improve convergence of the fixed-point iterations (which is often only linear and very slow) is to use acceleration methods [5,6]. In the specific case of the nucleation and growth problem, the concentration $A_5$ is the objective variable of the numerical methodology to solve the population balance equation, for this reason the acceleration method is applied only on the variable $A_5$. Indeed, the acceleration of a part of the output vector has demonstrated to be numerically efficient to converge strongly coupled multi-physics mathematical problems [7]. Three acceleration methods are tested to solve equation (2), namely: crossed secant, alternate secant and Steffensen. The last one is the only one to improve the convergence of the standard fixed-point iterations. In such a case, the accelerated iterate is defined as follows:
$$A_5^{j+1} = A_5^j - \frac{(f(A_5^j) - A_5^j)^2}{f(f(A_5^j)) - 2f(A_5^j) + A_5^j} \quad (3)$$

## 2.2 Convergence criteria

The population balance equation does not respond to a unique solution criterion. Indeed, several crystal size distributions $N_5$ can satisfy the same value of the limiting reagent concentration $A_5$. For this reason, the convergence criteria of the fixed-point algorithm is determined by the convergence in terms of concentration (equation (4)) and crystal size distribution (equation (5)).

$$|A_5^j - A_5^{j-1}| - \varepsilon_c < 0 \quad (4)$$

$$\max_{(k)} (|f(\mathbf{N}_5^j) - \mathbf{N}_5^j| - \varepsilon_a - \varepsilon_r |\mathbf{N}_5^j|) < 0 \quad (5)$$

with $\mathbf{N}_5^j$ the value of the iterate in the $j$th iteration, $A_5^j$ the concentration of the limiting reagent in the same iteration $\varepsilon_c$, $\varepsilon_a$ and $\varepsilon_r$ the convergence parameters. To converge, both equations (4) and (5) should be satisfied at the same iteration.

## 2.3 Initial condition

Iterative algorithms are particularly sensitive to the choice of a suitable initial condition. In this particular case, the iterating variable ($A_5$) is strongly coupled to the crystal size distribution ($N_5$) via the mass and population balances. For this reason, the initial condition concerns these two variables. The best initial condition is then given by the solution of both balances in the second reactor.

The recycle stream initial condition problem is then transferred to the initial condition of the inlet variables of the second reactor ($A_3$ and $N_3$ in the stream $F_3$ in Figure 1 b)). A previous work [3] demonstrated that the best initial condition to solve the population and mass balances in a homogenous reactor is a pure clear solution ($N_3^0 = 0$). Regarding the initial condition for the concentration, a zero value is not suitable because of the thermodynamic limit of the liquid-solid equilibrium. In addition, the initial condition cannot be equal to the solubility because this value implies a supersaturation equal to 1 and the absence of crystallisation phenomena. Consequently, the initial concentration ($A_3^0$) is chosen to be over the equilibrium value. Hence, the initial condition reads:

$$(A_5^0, N_5^0) = f_{R_2}(A_3^0, N_3^0), \quad N_3^0 = 0, \quad A_3^0 = \varepsilon_i * A_{eq}, \quad \varepsilon_i > 1 \quad (6)$$

Several values of the constant $\varepsilon_i$ were tested between 1.01 and 1.5. As a result, the convergence of the entire algorithm remains acceptable all over the interval. Also, $\varepsilon_i = 1.1$ represents the best final value of the concentration $A_3^0$.

# 3. Results

As a case study, the numerical methodology presented in the last section is applied to the uranium oxalic precipitation in a recycling structure):

$$U(NO_3)_4 + 2H_2C_2O_4 \xrightarrow{H_2O} U(C_2O_4)_2, 6H_2O + 4HNO_3 \quad (7)$$

At the industrial scale, the oxalic precipitation of uranium is carried out by adding simultaneously both reactants to the vortex reactor. However, in Figure 1, the stream $F_1$ represents the flowrate of both reactant streams. This simplification is merely schematic and does not modify the results of the modelling strategy.

The recycling structure is depicted in Figure 1b) and the volumes of the two reactors represent the hydrodynamics of the vortex reactor. Thus, the first reactor accounts for 5% of the total volume (forced vortex), the second is for 95% (free vortex) while the recycling ratio is 30%. Details on the crystallisation kinetics and the operating conditions can be found in the literature [3,8,9]. The numerical
methodology used to solve the mass and population balances is identical to the one developed for the precipitation of uranium oxalate [4].

## 3.1 Numerical analysis

Figure 2 presents the evolution of some state variables over the fixed-point iterations. By a simple inspection, the mass balance is satisfied first at the 5th iteration. Indeed, the variation of the uranium concentration $A_5$ is lower than $ε_c$, the imposed threshold. On the other hand, the population balance equation converges at about the 16th iteration even if the mean crystallite size remains stable from the 6th iteration. Indeed, the average size of the entire population is very close in two successive iterations but the entire crystallite size distribution does not match. The evolution of these three variables demonstrates the numerical bias imposed by the multiple solutions problem. The solution methodology described in the last section overcame this numerical issue by the choice of two simultaneous convergence criteria.

The numerical methodology developed in this paper is intended to be coupled with the agglomeration population balance equation and then applied to a more complex reactor network structure. In this concern, a low convergence time is required to contribute to the numerical performance of the entire algorithm. For the case presented in this paper, the convergence time is about 4 minutes in an Intel Core i7 machine (1.90 GHz/ 2.11 GHz) with 32 Go RAM.

<figure><img src="image_3.png" /><figcaption>Figure 2. Numerical analysis over the number of iterations: a) Concentration convergence criterion (equation (4)) b) Concentration of uranium c) Crystal size distribution convergence criterion (equation (5)) d) Mean crystallite size</figcaption></figure>

## 3.2 Process variables

The numerical methodology gives access to some valuable information on the recycle and the intermediate streams: Figure 3 depicts the crystal size distribution before and after the second reactor. For the same locations, the liquid phase composition, crystallisation kinetics, supersaturation, and suspension concentration are displayed in Table 1.

The mass transfer to the solid phase is mostly operated in the first reactor, less than 0.2 % of the final crystal mass is due to the crystallisation phenomena in the second reactor. Additionally, the average
crystal size increases by 0.2 µm and the shape of the crystal size distribution significantly changes. Indeed, the amount of small particles extremely decreases from 8*10<sup>23</sup> to 0 in the second reactor. This behaviour is explained by the relative values of the crystallisation kinetics in both homogeneous volumes. Nucleation and growth kinetics take lower values on the second reactor because of the drop in the supersaturation value. In fact, both heterogenous and homogeneous nucleation are negligible in the second volume, only the growth phenomenon takes place. Consequently, no crystal is born in the second reactor; the actinide precipitated is then consumed by growing the very small particles coming from the first reactor.

<figure><img src="image_4.png" /><figcaption>Figure 3. Crystallite size distribution.</figcaption></figure>

Another revealing information is the relationship between the residence time and the solid concentration. Most of the solid phase is created in a very small volume due to the nucleation process: 99.9% of the final uranium oxalate precipitates in 5% of the volume while only 0.1% occurs in 95% of the reactor volume. This demonstrates that the kinetics of the nucleation and growth processes are fast compared to the residence time. In addition, the volume of the reactor is large enough to bring the system until the thermodynamic limit imposed by the solubility of the uranium oxalate.

<table><caption>Table 1. Liquid and solid characteristics after algorithm convergence.</caption><thead><tr><th>Variables</th><th></th><th>R<sub>1</sub></th><th>R<sub>2</sub></th></tr></thead><tbody><tr><td>Residence time (s)</td><td>τ</td><td>50</td><td>955</td></tr><tr><td>Feeding ratio (-)</td><td>RS</td><td>0.36</td><td>0.006</td></tr><tr><td>Supersaturation (-)</td><td>S</td><td>6.3</td><td>1</td></tr><tr><td>Final concentration (mol m<sup>-3</sup>)</td><td>[U]</td><td>0.3</td><td>0.055</td></tr><tr><td></td><td>[C<sub>2</sub>H<sub>2</sub>O<sub>4</sub>]</td><td>53.9</td><td>53.3</td></tr><tr><td></td><td>[HNO<sub>3</sub>]</td><td>423.1</td><td>424.3</td></tr><tr><td>Solid concentration (mol m<sup>-3</sup>)</td><td>C<sub>c</sub></td><td>105.8</td><td>106</td></tr><tr><td>Growth kinetics (m s<sup>-1</sup>)</td><td>G</td><td>5.9 10<sup>-9</sup></td><td>1 10<sup>-12</sup></td></tr><tr><td>Homogeneous nucleation kinetics (m s<sup>-1</sup>)</td><td>R<sub>N</sub></td><td>9 10<sup>9</sup></td><td>0</td></tr><tr><td>Heterogeneous nucleation kinetics (m s<sup>-1</sup>)</td><td>R<sub>N</sub></td><td>4.8 10<sup>15</sup></td><td>0</td></tr><tr><td>Crystals average size (µm)</td><td>d<sub>4,3</sub></td><td>1.6 10<sup>-6</sup></td><td>1.8 10<sup>-6</sup></td></tr></tbody></table>
# 4. Conclusion and perspectives

The study presented in this paper shows that the accelerated fixed-point algorithm appears to be adapted to the solution of the recycle stream in the reactor network approach. The resulting numerical methodology predicts the crystal size distribution and actinide concentration in the context of the uranium oxalate precipitation.

The solution of the recycle stream coupled with the population balance solution methodology gives valuable information about the hierarchy between the crystallisation phenomena at the local scale. At industrial scale, this information can be used to simultaneously adjust the crystal size distribution and the actinide concentration in the product stream.

The numerical methodology introduced in this study is easily adaptable to deal with any distributed variable over a recycle stream, independently of the numerical strategy adopted to solve the equations system in the homogenous volume. Future work will include the solution of the agglomeration phenomenon over the recycle stream.

# Acknowledgments

The authors thank Orano for its financial support.

# References

[1] N. Jourdan, T. Neveux, O. Potier, M. Kanniche, J. Wicks, I. Nopens, U. Rehman, and Y. Le Moullec, Chem. Eng. Sci. 210, 115196 (2019)

[2] E. Saikali, M. G. Rodio, G. Bois, U. Bieder, N. Leterrier, M. Bertrand, and Y. Dolias, Chem. Eng. Sci. 214, 115426 (2020)

[3] C. C. Ruiz Vasquez, N. Lebaz, I. Ramière, S. Lalleman, D. Mangin, and M. Bertrand, Chem. Eng. Res. Des. 177, 767 (2021)

[4] C. C. Ruiz Vasquez, N. Lebaz, I. Ramière, S. Lalleman, D. Mangin, and M. Bertrand, J. Cryst. Growth 616, 127258 (2023)

[5] I. Ramière and T. Helfer, Comput. Math. Appl. 70, 2210 (2015)

[6] C. Brezinski and L. Wuytack, Numerical Analysis: Historical Developments in the 20th Century. (Elsevier Science, Amsterdam, 2014)

[7] Isabelle, Ramière, Autour de Méthodes Numériques Pour Les Couplages Multiphysiques et Multiéchelles En Mécanique Des Solides., Aix Marseille Université (AMU), 2021

[8] M. Bertrand-Andrieu, E. Plasari, and P. Baron, Can. J. Chem. Eng. 82, 930 (2008)

[9] S. Lalleman, M. Bertrand, and E. Plasari, J. Cryst. Growth 342, 42 (2012)