Application of Computational Fluid Dynamics for Modeling of Hydrodynamics and Mass Transfer of Laboratory Scale Crude Palm Oil Degumming Process

In this research, computational fluid dynamics was applied to model a laboratory reactor for the degumming of crude palm oil (CPO) involving immiscible liquid-liquid mixing and being controlled by mass transfer. The fluid mixing of CPO and phosphoric acid in the reactor was simulated using multiphase mixture model and standard k-ε turbulence model in steady state mode. The simulation predicts the distributions of the drop diameter, the dispersed-phase relative velocity, the drop Reynolds number and the interfacial area density. The mass transfer coefficient from experimental work is correlated using the model as Sherwood number Shd = 0.02576 RedSc with R being 0.91. Finally, the volumetric mass transfer coefficient was calculated and compared to the experiment result.


Introduction
Degumming is an important step in crude palm oil (CPO) refining. It reduces or removes gum content from CPO. Gum consists of phospholipid compounds that acts as an emulsifier. The existence of gum disrupts product separation [1].
One of degumming methods, called acid degumming, uses phosphoric acid solution (H 3 PO 4 ) as a reagent [2]. Non-hydratable phospholipids in CPO reacts with the solution and changes into hydratable phospholipids which is easy to remove [3].
The mechanism of CPO degumming involves the phospholipid transfer from the CPO bulk to the H 3 PO 4 solution through the interface, followed by the chemical reaction of phospholipids and H 3 PO 4 . The overall rate of degumming depends on the mass transfer coefficient and the reaction rate coefficient. The volumetric mass transfer coefficient is a function of the total mass transfer coefficient and the interface area which are usually promoted by agitation in a reactor. Therefore, the mixing effectiveness is highly influential to the mass transfer rate.
Several methods have been proposed to correlate mixing dynamics in a reactor to mass transfer rates. Some correlations are merely based on dimensional analysis and thus fully empirical, while others are based on turbulence and relative velocity theories. The latter is built upon the effect of forced convection around particles on mass transfer rate. Other authors derived expressions for relative velocities by solving motion equations of suspended particles in flow fields [4]. This approach uses local parameters (relative velocities between phases) to describe flow dynamicsdependentmass transfers. The approach evolves the differences in flow fields in different geometries of reactors.
Relative velocities can be calculated using semianalytical methods with some assumptions. Therefore, it is limited by flow geometries of mixing process. Otherwise, it can be determined by solving multiphase flow equations using numerical techniques, which is called computational fluid dynamics (CFD).
CFD uses differential fluid flow governing equations and is solved numerically. It predicts flow parameters of a system, regardless of its geometry. Due to its flexibility, CFD has been widely applied for modeling and designing many kinds of reactors and processes. For liquid-liquid systems, CFD is used by applying multiphase Eulerian models and population balances to predict liquid-liquid mixing and mass transfer in a reactor [5,6]. A simple mixture model has been used as substitute for more complicated Eulerian model to design a submerged membrane stirred reactor [7].
In multiphase mixture models, continuity equations are solved for each phase. On the other hand, momentum equations are solved for the mixture of continuous and dispersed phases in case of two-phase flows. The model is valid where the local equilibrium between the phases is assumed. When dispersed phase clustering and large slip velocity are possible, the model would be invalid [8]. The two conditions usually arise in mixing with large dispersed phase volume fraction (> 20%).
In the present study, a mixture model is used to simulate crude palm oil (CPO) degumming process in a laboratory reactor. The simulation result, i.e. the drop diameter distribution, was compared to empirical data to obtain the interfacial area. Mass transfer coefficient was then calculated and correlated to mixing flow parameters. The model was also solved to predict the flow distribution and pattern so that the effect of different internal geometries of the reactor can be quantified.

Laboratory reactors
The laboratory data used in the study is obtained from Ristianingsih et al. [9]. The volume of the reactant mixture was 100 mL consisting of CPO as continuous phase and phosphoric acid solution (H 3 PO 4 85%) as dispersed phase. The volume fraction of dispersed phase was 1.5%. The reaction took place in a 1000 mL tripleneck round vessel. A magnetic stirrer bar was used to promote agitation. The temperature was kept constant during the reaction. The reaction was run for 120 minutes with initial content of gum in CPO being 1.43%. The parameters deduced from the measurement were the volumetric mass transfer coefficient (k c a), the reaction rate coefficient (k) and the equilibrium constant (K). The set of experiment data is shown in Table 1, where k c a becoming the concern. The properties of each component used in the simulation is summarized in Table 2. For this study, all properties were not measured. They were calculated using empirical equation available in literatures.

Geometrical modeling
The modeling requires the drawing of the system geometry. The fluid domain of the laboratory reactor was drawn based on its internal configuration. The reactor internal contains only a magnetic stirrer bar submerged in the reactant fluid, as depicted in Figure 1. The fluid is divided into two domains: static and rotating. The rotating domain contains the movable reactor internal.

Mathematical modeling
The mixing inside the reactor is predicted by solving the differential equations: continuity, momentum balance, mixture model, slip velocity and turbulence k-ε. In order to model the rotating fluid inside the reactor, the multiple reference frames (MRF) equation is applied.
In the mixture model, the continuity and momentum equations are in the form as follow [10]: The velocities of each phase are calculated by using the slip velocity equation. In this case, the Hadamard-Rybczynzki model is used: The drag coefficient is calculated as follow: The turbulence model is used to predict the fluid flow instability generated by the stirrer bar rotation. It uses the standard k-ε model. The turbulence is manifested as additional stress in the momentum equation. This stress is calculated by solving two transport equations, i.e. the turbulence kinetic energy and the turbulent dissipation rate: The rotating motion inside the reactor is described using MRF equation by introducing additional forces to the momentum equation, i.e. coriolis and centrifugal forces.
In order to evaluate the distributions of the dispersed-phase drop diameter and the interface area, the transport equations for the dispersed-phase volume fraction and the drop number density need to be solved. The transport of the two quantities is determined by the velocity of the dispersed phase obtained from Equation (1) to Equation (6). After several algebraic modifications, the two transport equations become The diameter and the interface area then can be calculated using following equations: The two equations take the assumption that the number of the dispersed phase drops is constant throughout the mixing process.
The drop diameter calculated using Equation (9) was compared to that calculated using an empirical equation derived from direct measurement where the data was gathered from published work [11]. The method of developing measurement data-basedempirical equations has been presented by Rueger and Calabrese [12]. The data is better to be classified based on flow regimes, because the dispersion mechanism is different. The effect of dispersed-phase volume fraction has to be also considered.
Based on the plotting of impeller power numbers versus impeller Reynolds numbers, the mixing flow is discovered within laminar regime to intermediate regime. The forces acting on the dispersed phase drops are predominantly in the form of shear stress. Therefore, the empirical equation to estimate the average diameter (d 32 ) of the drops is C a is capillary number as the function of viscosity ratio (λ), which can be derived from the plot published by Grace [13].
Equation (11) is applied for a dilute mixture with dispersed-phase volume fraction being ≤ 0.1%. To adjust the estimate of average drop diameter due to the effect of volume fraction, the correction factor has been introduced by Rueger and Calabrese [14].

CFD simulation
The model wassolved using COMSOL Multiphysics. The physics used was Rotating-Machinery Mixture-Model Turbulent Flow. Simulation run in stationary mode (Frozen Rotor). This mode was used to simulate the mixing at steady state. Simulation was conducted for all operating condition mentioned in Table 1.

Results and discussion
Several variables solved by simulation are velocity of mixture, pressure, slip velocity, dispersed phase volume fraction, and number density. From this solution, other variables can also be derived, include drop diameter, interfacial area, drop relative velocity, and drop Reynolds number.  From the streamline pattern, it seems that solid body vortices are formed following the direction of stirrer bar rotation. The smaller vortex is also generated near the wall. The fluid flow in vertical section, referring to Fig.   2(b), directs toward the reactor wall as the effect ofa centrifugal force and then returns to the center of stirrer bar rotation. A similar flow pattern was also pointed out by previous study, supported by LDV (laser Doppler velocimetry) measurement [15].

Mixing hydrodynamics
The turbulence is not similarly perceived by the drops throughout the reactor. Several authors proposed the usage of drop relative velocityto define local turbulence [4]. This is manifested as drop Reynolds number defined as: The drop relative velocity andthe drop diameter are available as simulation results. Therefore, the distribution of drop Reynolds number can be calculated.
It is showed in Fig. 3giving the picture of local turbulence inside the reactor.

Fig. 3. Distribution of drop Reynolds number for case 500 rpm, 80 °C
The turbulence is more intense around the stirrer bar tip than near the reactor wall. This might be associated with the momentum generated by stirrer bar motion. This momentum is then dissipated away from the rotation zone, also with the turbulence intensity.

Mass transfer
The purpose of this modeling is to generate a correlation between the local turbulence intensity to the mass transfer coefficient. This correlation, then coupled with interfacial area predicted from simulation, is used to estimate the distribution of the volumetric mass transfer coefficient (k c a). However, due to the availability of the laboratory reactor data as a single value instead of a distribution, the drop Reynolds number value has to be averaged. The averaging method used is the geometrical averaging. This is done by considering the statistical distribution of the variable that approaches the lognormal distribution. Moreover, due to several laboratory experimentsis also performed in different temperature, its effect is also involved in the correlation.
To generate a correlation, each of above factors is represented as dimensionless number, i.e. Sherwood number, drop Reynolds number, and Schmidt number. By applying multivariable regression, an equation to correlate the dimensionless numbers is obtained: The coefficient of determination (R 2 ) value of 0.91 shows good correlation between the model and the data asseen in Fig. 4.

Fig. 4. Comparison between Sherwood number from data and correlation
From Equation (13), due to its higher exponent,the drop Reynolds number has stronger influence to the Sherwood number rather than to Schmidt number. Therefore, the mass transfer coefficienthas more dependence to mixing intensityrather than to temperature.
The second factor affecting the volumetric mass transfer coefficient is interfacial area. This variable is derived from the simulation result by Equation (10). The distribution of the interfacial area density is presented in state. The dropsare gathered in the front of the stirrer bar tip during the rotation so that the number density and the volume fraction are higher. This leads tohigher interfacial area density. Behind the stirrer bar tip, the effect reverses. This happensbecause the turbulence and therefore the dispersion are more intensive by the existence of trailing vortex [16].  The k c a is higher near the stirrer bar rotation zone, especially near its tip, and is lower approaches the near wall zone. Therefore, the effective zone for the mass transfer to occur is located around the stirrer bar rotation. Fig. 6 (b) shows the statistical distribution of k c aderived from the simulation. Its average value is compared to the k c a obtained from the laboratory experiments. The relative error of both values is around 20%. This error is probably related to the error in the Sherwood number correlation. Although the error is rather high, from the modeling point of view, the estimation is quite reasonable.This is because it combines the fundamental equation of fluid flow and the experimental data.Using this modeling approach, value of k c a for different reactor configurations (e.g. with different impellers or reactor sizes) can be estimated prior to actual experiments. This is useful, for example, if a reactor scale-up is planned.
The CFD simulation has been carried out for the laboratory reactor of CPO degumming process. The simulation, combined with the data analysis of the laboratory experiments and the empirical equation to estimate the drop size, is used to generate Sherwood number correlation. The correlation shows the good coefficient of determination (R 2 ) of 0.91. By applying this correlation and coupled with the value of the interfacial area from the simulation, the distribution of the volumetric mass transfer coefficient (k c a) in the reactor is resulted. The average of this k c a is comparable to k c a from laboratory experiments with the relative error being ~20%. This modeling approach is useful in estimating k c a of other CPO degumming reactor with different sizes or configurations.