CFD Modeling of Spatial Inhomogeneities in a Vegetable Oil Carbonation Reactor

: Fossil materials are widely used raw materials in polymerization processes; hence, in many cases, the primary goal of green and sustainable technologies is to replace them with renewables. An exciting and promising technology from this aspect is the isocyanate-free polyurethane production using vegetable oil as a raw material. Functional compounds can be formed by the epoxidation of vegetable oils in three reaction steps: epoxidation, carbonation, and aminolysis. In the case of vegetable oil carbonation, the material properties vary strongly, with the composition a ﬀ ecting the solubility of CO 2 in the reaction mixture. Many attempts have been made to model these interactions, but they generally do not account for the changes in the material properties in terms of spatial coordinates. A 2D CFD model based on the combination of the k- ε turbulence model and component mass balances considering the spatial inhomogeneities on the performance of the reactor was created. After the evaluation of the mesh independence study, the simulator was used to calculate the carbonation reaction in a transient analysis with spatial coordinate-dependent density and viscosity changes. The model parameters (height-dependent mass transfer parameters and boundary ﬂux parameters) were identiﬁed based on one physical experiment, and a set of 15 experiments were used for model validation. With the validated model, the optimal operating temperature, pressure, and catalyst concentration was proposed. case of experiment it is 413.15 K. The results are shown at the end of the simulation time. As we can see, the velocity ﬁeld of experiments 8 and 12 are di ﬀ erent, meaning that the temperature has a more signiﬁcant e ﬀ ect on the velocity ﬁeld than pressure.


Introduction
Fossil raw materials are the most commonly used raw materials in the energy sector and the chemical industry. Nearly 90% of organic matter-based products are produced from petroleum or natural gas. The cornerstone of sustainable technology development is to reduce the dependence on fossil raw materials, which is an essential directive at the industrial and governmental levels as well. Using vegetable oil as a renewable raw material is an excellent example of this effort. The most important characteristics of vegetable oil are its accessibility and excellent biodegradability. Different feedstocks can be used for the process-for instance, soybean, rapeseed, sunflower seed, palm oil etc. [1]. The properties of the vegetable oil vary with the length of the chains and the number of double bonds. The application areas of vegetable oils are vast; they are used to produce different type of polyurethanes, resins, and even composites, which make vegetable oils promising compounds in the future of chemistry [2]. One of the most promising fields can be the production of secondary fuel for diesel vehicles [3,4].
In this study, we focus on the production of widely used isocyanate-free polyurethane (this polymer is the sixth most commonly used material). Several other studies dealt with the investigation of 2.1. Geometry Selection, Momentum Balance Calculation for the 2D and the 3D Model The reactor was a Parr 300 mL reactor with 100 mm inner diameter equipped with a gas entrainment impeller (25 mm diameter). At the beginning of each experiment, 100 mL of epoxidized cottonseed oil (ECSO) was loaded into the reactor. Further information about the measurements can be found in Cai et al. [5]. A 500 rpm revolution speed was applied in each experiment. A 3D and a 2D representation of the reactor geometry were implemented in the CFD model. In the case of the 3D representation, a rotating reference frame method was used, while swirl flow was implemented in the case of the 2D representation. In this study, an axisymmetric 2D representation of the reactor was implemented in COMSOL Multiphysics. The main idea of the axisymmetric calculation is the facilitation of the symmetry boundaries; thus, we can decrease the computational time significantly. In this case, the reactor itself is not symmetric due to the asymmetry of the impeller. However, the relatively fast revolution speed (500 rpm) makes the flow field relatively symmetric, which justifies using this approach. Figure 1a shows the 3D and the 2D representation, and Figure 1b shows the resulting velocity fields (m/s). The velocity fields were similar in the case of 3D and the 2D simulations. However, a multiplication factor of 2.8 has to be implemented in the revolution speed in case of the 2D simulation. To make sure that the numerical values are similar in the case of the axial symmetry model, a reasonable agreement was found between the velocity field of the 3D and the 2D representation, so the 2D model was used for the further parameter identification and optimization steps. The main advantage of using the 2D model is the lower computation time, making the ordinary optimization methods applicable.  Figure 1a shows the 3D and the 2D representation, and Figure 1b shows the resulting velocity fields (m/s). The velocity fields were similar in the case of 3D and the 2D simulations. However, a multiplication factor of 2.8 has to be implemented in the revolution speed in case of the 2D simulation. To make sure that the numerical values are similar in the case of the axial symmetry model, a reasonable agreement was found between the velocity field of the 3D and the 2D representation, so the 2D model was used for the further parameter identification and optimization steps. The main advantage of using the 2D model is the lower computation time, making the ordinary optimization methods applicable.  Figure 2a shows the reactants (epoxidized cottonseed oil, ECSO) and the products (carbonated cottonseed oil, CCSO). Figure 2b presents the reaction mechanism. In this work, the model of Cai et al. was extended to consider spatial inhomogeneities. Later on, we will refer to this model as the original one [5].  Figure 2a shows the reactants (epoxidized cottonseed oil, ECSO) and the products (carbonated cottonseed oil, CCSO). Figure 2b presents the reaction mechanism. In this work, the model of Cai et al. was extended to consider spatial inhomogeneities. Later on, we will refer to this model as the original one [5].

Carbonation Reaction Modeling
The viscosity of the reaction mixture, which can be calculated based on the sum of weighted viscosities of the components varies in the spatial and temporal coordinates. Equation (1) presents how the reaction rate is calculated. The viscosity of the reaction mixture, which can be calculated based on the sum of weighted viscosities of the components varies in the spatial and temporal coordinates. Equation (1) Figure 2. The applied reaction mechanism (a) simplified (b) detailed with catalyst [5].
The maximum solubility of CO 2 (c CO2_liq * ) is calculated within the CFD simulator based on the Henry coefficients introduced in the original model. Equations (2) and (3) show the calculation of c CO2_liq* and the overall He coefficient.
He mixing = x ECSO · He ECSO + x CCSO · He CCSO Equation (4) shows the first addition to the original model, the flux of the CO 2 at the inlet boundary, which was the contact area between the fluid and the gas phase. The component viscosities and densities are calculated with the following Equations (5)-(8), based on the measurements described in Cai et al. [5].
The material properties were implemented using Equations (9) and (10), including temperature and coordinate-based calculation. Another extension of the original model is the calculation of the N CO2 term, where the mass transfer coefficient and the vertical change in the bubble size were considered, where z refers to the spatial z coordinate. A and B are the constants to be identified.
Equations (12)-(14) present the component balances, which are in similar forms as in the original model.

CFD Model of the Reactor
One of the main goals was to calculate the spatial and temporal changes of the viscosity and the density. In most of the CFD simulators, when the momentum balance is close to stationary (or reaches the stationary state significantly faster than the component balance) the two-equation sets can be separated, and a stationary momentum balance can be used as a basis for the calculation of the component balances. In our case, the viscosity and the density, which are the two main material properties used in the momentum balance, change in time. Since they vary with the progression of the reaction, the two balances cannot be separated. Thus, a longer computational time is expected.
The reactor was operated isothermally, so we neglected the solution of the heat balances. A mesh independence study was performed to ensure the validity of the model. Five different mesh sizes were applied, and the changes in the velocity field were calculated within a grid. The changes in the velocity field are getting lower by the increase of the number of mesh elements, and consequently, it results in a higher computational time. After the 4th mesh size, the changes in the velocity become around 1%, while the computational time tripled from 4th to 5th mesh, so the 4th mesh size was chosen for the further calculations. Figure 3 shows the results of the mesh independence study and the resulting mesh.
The reactor was operated isothermally, so we neglected the solution of the heat balances. A mesh independence study was performed to ensure the validity of the model. Five different mesh sizes were applied, and the changes in the velocity field were calculated within a grid. The changes in the velocity field are getting lower by the increase of the number of mesh elements, and consequently, it results in a higher computational time. After the 4th mesh size, the changes in the velocity become around 1%, while the computational time tripled from 4th to 5th mesh, so the 4th mesh size was chosen for the further calculations. Figure 3 shows the results of the mesh independence study and the resulting mesh.

Experimental Work Considered
Based on a previous article [5], 14 experiments were considered for this study. The investigated temperature interval was between 110 and 140 °C; the pressure was between 30.6 and 52 bar. Table 1 presents the operating conditions of the experiments.

Experimental Work Considered
Based on a previous article [5], 14 experiments were considered for this study. The investigated temperature interval was between 110 and 140 • C; the pressure was between 30.6 and 52 bar. Table 1 presents the operating conditions of the experiments.

Results
This section introduces the results of the model parameter identification (see Section 3.1) and the results of the optimization problem (see Section 3.3). Section 3.2 describes the sensitivity analysis of the revolution speed based on the first experiment.

Model Parameter Identification
Three parameters were identified, which are the A and B height-dependent mass transport parameters in Equation (12), and the D component transfer parameter from Equation (4) using experiment 1 (see Table 1). The NOMAD black-box optimization toolbox was used for parameter identification [8,9]. The following objective function (Equation (15)) is applied to identify the values of the A, B, and D parameters, which is a sum of the relative squared errors in time.
where k is the number of measurement points (in time).
The upper and lower limits for the parameters A, B, and D were (3, 15, 3 × 10 −4 ) and (1, 10, 10 −4 ), respectively. A COMSOL-MATLAB Livelink connection was used for the calculation, where the 2D CFD model was calculated in every function call. An Intel Xeon E5620 computer was used for the calculation with 72 GB RAM. One function call lasted around five minutes. A grid size-based integrated concentration was used to compute the simulated value. The parameter identification step resulted in A = 2.101, B = 13.223, and D = 1.984 × 10 −4 (m/s). The summed error was 9.45. Figure 4 shows the integrated ECSO concentrations. The markers show the measured values for each experiment, while the continuous lines show the simulated values.
As can be seen in Figure 4, most of the simulated values are in good agreement with the measured values. To further evaluate the effect of the process parameters of the errors, the connection between the temperatures, pressures, and process parameters are plotted in Figure 5.
As we can see, the model error is higher at lower temperatures (see Figure 4; exp. 2 and 12). The higher pressures correlate with higher errors (see Figure 4; exp. 5, 11, 12 and 13).
In the next section, the effect of the temperatures and pressures are shown on the velocity field and the concentrations profiles. Firstly, we show the results of experiments 3 and 13, which have the same temperature, and then experiments 8 and 12, which have similar pressure. Figure 6 presents the velocity field and the concentration of ECSO and CO 2 in the case of experiments 3 and 13. The temperature is the same in both experiments, but the pressure of experiment 13 is twice as much as in experiment 3. The velocity fields and ECSO and CO 2 concentrations are shown at the end of the simulation time.
The velocity field is different only in the maximum velocity; the higher pressure will lead to higher maximum velocity (0.721 compared to 0.717) as well as higher CO 2 concentration. The concentrations fields are relatively homogeneous.   As we can see, the model error is higher at lower temperatures (see Figure 4; exp 2 and 12). The higher pressures correlate with higher errors (see Figure 4; exp 5,11,12 and 13).    As we can see, the model error is higher at lower temperatures (see Figure 4; exp 2 and 12). The higher pressures correlate with higher errors (see Figure 4; exp 5,11,12 and 13).  experiments 3 and 13. The temperature is the same in both experiments, but the pressure of experiment 13 is twice as much as in experiment 3. The velocity fields and ECSO and CO2 concentrations are shown at the end of the simulation time. The velocity field is different only in the maximum velocity; the higher pressure will lead to higher maximum velocity (0.721 compared to 0.717) as well as higher CO2 concentration. The concentrations fields are relatively homogeneous.   Figure 7 shows the velocity field and the concentrations of ECSO and CO 2 in the case of experiments 8 and 12. The pressure is similar, but the temperature of experiment 12 is 383.15 K, while in the case of experiment 8, it is 413.15 K. The results are shown at the end of the simulation time. As we can see, the velocity field of experiments 8 and 12 are different, meaning that the temperature has a more significant effect on the velocity field than pressure. experiments 8 and 12. The pressure is similar, but the temperature of experiment 12 is 383.15 K, while in the case of experiment 8, it is 413.15 K. The results are shown at the end of the simulation time. As we can see, the velocity field of experiments 8 and 12 are different, meaning that the temperature has a more significant effect on the velocity field than pressure.
The ECSO and the CO2 concentration profiles get more homogenous with the increase of temperature compared to the lower temperature cases. The viscosity and the densities are time and spatial coordinate dependent. Figure 8 shows the viscosity and density changes in the function of time in the case of experiment 1. As we can see, the dynamic viscosity of the reaction phase changes by around 7%, while the density change is around 3%. The changes in time are significant, but the changes in space are minimal due to adequate mixing. The ECSO and the CO 2 concentration profiles get more homogenous with the increase of temperature compared to the lower temperature cases.
The viscosity and the densities are time and spatial coordinate dependent. Figure 8 shows the viscosity and density changes in the function of time in the case of experiment 1. As we can see, the dynamic viscosity of the reaction phase changes by around 7%, while the density change is around 3%. The changes in time are significant, but the changes in space are minimal due to adequate mixing.

Sensitivity Analysis of the Revolution Speed (Experiment 1)
In this section, we examined the effect of the different revolution speeds. The calculation was completed with the simulator, including the component mass balances.
In addition to the initial case, two different revolution speeds were applied (100 rpm and 300 rpm). Figure 9a-c shows the velocity fields.

Sensitivity Analysis of the Revolution Speed (Experiment 1)
In this section, we examined the effect of the different revolution speeds. The calculation was completed with the simulator, including the component mass balances.
In addition to the initial case, two different revolution speeds were applied (100 1/min and 300 1/min). Figure 9a-c shows the velocity fields. As we can see, not only the magnitude but the characteristics of the flow are changing with the increase of the revolution speed. The changes of ECSO concentration were minimal. The change of the velocity magnitude is not linear dependent on the revolution speed; however, to extend the model for operating with different revolution speeds, additional measurements are needed.

Optimization
The reactor performance was optimized, considering the three most critical process parameters: the process temperature, the pressure, and the initial catalyst concentration. The NOMAD black-box optimization algorithm was applied to solve the optimization problem. COMSOL-MATLAB Livelink was used for the solution. The CFD simulator was implemented in COMSOL, while the optimization algorithm was implemented in a MATLAB environment. The upper and the lower limits of the optimization problem were (413.15 60 0.7) and (383.15 30 0.5), respectively. The temperature and pressure boundaries were defined based on the physical limitations of the systems. The objective is to maximize the concentration of the carbonated cottonseed oil (CCSO) at the end of the reaction. The optimal parameters were found as 383.39 K and 51.14 bar, while the optimal initial concentration of the catalyst is 0.699 mol/L. The outlet concentration of ECSO is 1.36 × 10 −5 mol/L. Figure 10 shows the integrated concentrations of epoxidized cottonseed oil (ECSO) and the carbonated cottonseed oil (CCSO). As we can see, not only the magnitude but the characteristics of the flow are changing with the increase of the revolution speed. The changes of ECSO concentration were minimal. The change of the velocity magnitude is not linear dependent on the revolution speed; however, to extend the model for operating with different revolution speeds, additional measurements are needed.

Optimization
The reactor performance was optimized, considering the three most critical process parameters: the process temperature, the pressure, and the initial catalyst concentration. The NOMAD black-box optimization algorithm was applied to solve the optimization problem. COMSOL-MATLAB Livelink was used for the solution. The CFD simulator was implemented in COMSOL, while the optimization algorithm was implemented in a MATLAB environment. The upper and the lower limits of the optimization problem were (413.15 60 0.7) and (383.15 30 0.5), respectively. The temperature and pressure boundaries were defined based on the physical limitations of the systems. The objective is to maximize the concentration of the carbonated cottonseed oil (CCSO) at the end of the reaction. The optimal parameters were found as 383.39 K and 51.14 bar, while the optimal initial concentration of the catalyst is 0.699 mol/L. The outlet concentration of ECSO is 1.36 × 10 −5 mol/L. Figure 10 shows the integrated concentrations of epoxidized cottonseed oil (ECSO) and the carbonated cottonseed oil (CCSO). Processes 2020, 8, x FOR PEER REVIEW 13 of 16 Figure 10. Integrated concentrations of ECSO and CCSO with the optimal operating parameters 383.39 K, 51.14 bar, 0.699 mol/L.
As we can see, full conversion can be achieved with the optimal parameters. Figure 11 shows the velocity field and the concentration of ECSO within the reactor after 4 h. The velocity field is similar in low-temperature cases. The ECSO concentration is reduced to zero almost in the whole reactor volume, only a minimal amount remains. However, the main driving force for the higher conversion is the higher concentration of CO2 in the mixture, which is the result of the higher pressure. Figure 12 shows the changes in viscosity and density in time. As the reaction progresses further than the experimental cases, the changes of viscosity and density are higher as well. As we can see, full conversion can be achieved with the optimal parameters. Figure 11 shows the velocity field and the concentration of ECSO within the reactor after 4 h. As we can see, full conversion can be achieved with the optimal parameters. Figure 11 shows the velocity field and the concentration of ECSO within the reactor after 4 h. The velocity field is similar in low-temperature cases. The ECSO concentration is reduced to zero almost in the whole reactor volume, only a minimal amount remains. However, the main driving force for the higher conversion is the higher concentration of CO2 in the mixture, which is the result of the higher pressure. Figure 12 shows the changes in viscosity and density in time. As the reaction progresses further than the experimental cases, the changes of viscosity and density are higher as well. The velocity field is similar in low-temperature cases. The ECSO concentration is reduced to zero almost in the whole reactor volume, only a minimal amount remains. However, the main driving force for the higher conversion is the higher concentration of CO 2 in the mixture, which is the result of the higher pressure. Figure 12 shows the changes in viscosity and density in time. As the reaction progresses further than the experimental cases, the changes of viscosity and density are higher as well.   Figure 12 shows the spatial coordinate-dependent changes in density and viscosity in case of the optimum. The viscosity increase is 35%, while the density increase is around 5%. As we can see, the viscosity depends on the ratio of ECSO/CCSO. In the case of lower temperatures and high pressure, the viscosity changes affect mainly the velocity field. The developed methods can serve well even with lower revolution speeds; however, to validate low revolution speed results, additional measurements will be needed.

Conclusions
In this study, a 2D CFD model of a carbonation reactor was developed. The model was based on an axial symmetric approximation of the system using the swirl flow model after the comparison of the 3D counterpart. The turbulent momentum balance was calculated together with the component mass balances in a transient simulation calculating the spatial coordinate-dependent changes of viscosity and density in the function of the component concentration, which is the main contribution of this study to the field.
The model parameters were identified based on one experiment and validated against 13 distinct measurements minimizing a relative squared error between the experimental and simulation results. The changes in concentrations, velocity fields, and the material parameters are discussed.
Then, the validated simulator is applied to optimize the optimal operating parameters, which are 383.39 K, 51.14 bar, and 0.699 mol/L concentration of the catalyst.