A Theoretical Analysis on a Multi-Bed Pervaporation Membrane Reactor during Levulinic Acid Esterification Using the Computational Fluid Dynamic Method

Pervaporation is a peculiar membrane separation process, which is considered for integration with a variety of reactions in promising new applications. Pervaporation membrane reactors have some specific uses in sustainable chemistry, such as the esterification processes. This theoretical study based on the computational fluid dynamics method aims to evaluate the performance of a multi-bed pervaporation membrane reactor (including poly (vinyl alcohol) membrane) to produce ethyl levulinate as a significant fuel additive, coming from the esterification of levulinic acid. For comparison, an equivalent multi-bed traditional reactor is also studied at the same operating conditions of the aforementioned pervaporation membrane reactor. A computational fluid dynamics model was developed and validated by experimental literature data. The effects of reaction temperature, catalyst loading, feed molar ratio, and feed flow rate on the reactor’s performance in terms of levulinic acid conversion and water removal were hence studied. The simulations indicated that the multi-bed pervaporation membrane reactor results to be the best solution over the multi-bed traditional reactor, presenting the best simulation results at 343 K, 2 bar, catalyst loading 8.6 g, feed flow rate 7 mm3/s, and feed molar ratio 3 with levulinic acid conversion equal to 95.3% and 91.1% water removal.


Introduction
Global warming caused by the growing greenhouse gases generation has been currently recognized as a major environmental problem. The use of petroleum-derived fossil fuel and the concomitant vegetation areas reduction are the main causes of a rise in CO 2 concentration in the atmosphere. Apart from the CO 2 effect, an increase in other industrial byproducts such as chlorofluorocarbons, hydrocarbons, methane, nitrogen oxides (NOx), and sulfur oxides (SOx) also increase the greenhouse effect. From an environmental point of view, the main advantage of using biofuels is reducing the effect of greenhouse and acid rain. As a renewable fuel option, in the near future, the importance of using biodiesel and bio-ethanol and, in the medium future, the importance of using biomass-derived Fischer-Tropsch diesel and cellulosic biomass-based additives will increase. Biomass-based ethyl levulinate (EtLA) is a significant fuel additive. The use of EtLA as a bio-fuel provides higher engine efficiency and lower CO and NOx emissions in addition to providing long operating life. Since diesel engines are responsible for high NOx and exhaust gas emissions, they have negative effects on the environment. To solve this issue, one of the most effective methods consists of the addition of oxygenates to fossil fuels, which are known as alkyl levulinate as a fuel oxygenate. EtLA attracted special attention due to its properties of high ignition temperature, 33% oxygen content, and high-efficiency clean combustion [1][2][3][4][5][6]. Levulinate esters are important chemical compounds used as fuel additives, solvents, and plasticizers. Ethyl levulinate can be used as a bio-additive in diesel fuel up to 5 wt%. The EtLA is synthesized by esterification of levulinic acid (LA) with ethanol (Et) in the presence of homogeneous or heterogeneous acid catalysts. Both reactants are produced from biobased sources. The use of bio-based reactants makes the production of EtLA process green and environmentally friendly [7,8]. EtLA is usually synthesized by using acid catalysts such as H 2 SO 4 , HCl, and H 3 PO 4 . They are unrecyclable and corrosive. Therefore, in recent years, many research groups explored the potential of green catalysts use. Heterogeneous solid acid catalysts are a suitable alternative to overcome the drawbacks of homogeneous acid catalysts. Heterogeneous catalysts can be separated from the reaction mixtures easily and reused for repeated runs [8][9][10][11].
Nowadays, membrane processes belong to a large group of techniques for separating the components of liquid and gas mixtures, depending on the properties of the used membranes and the specific characteristics of the separation itself. Different membrane operations are currently adopted for separating products ranging in size from tens of µm to tenths of nm, and particular types of membranes allow for the separation of vapors, gases, ions, etc. (reverse osmosis, nanofiltration, and membrane distillation). Among them, membrane pervaporation (PV) is useful for separating liquid mixtures, including those forming azeotropes. It depends on the selective sorption and diffusion rate differences and not from the liquid-vapor equilibrium. Most of the polymeric membranes used currently in PV processes may be categorized as in the following: • Inert membranes possessing hydrophilic nature (i.e., polyvinyl alcohol, cellulose acetate); • Inert membranes possessing hydrophobic nature (i.e., polydimethylsiloxane); • Hydrophilic ion exchange membranes.
Furthermore, membrane engineering deals also with the combination of two or more operation units in order to provide benefits such as the reduction in the cost of investment and the energy consumption. In this regard, pervaporation membrane reactors (PVMRs) may perform both the reaction and separation stages simultaneously. The application of PVMR enhances the conversion of reversible reactions as a result of the continuous removal from the reaction mixture of the byproduct with a membrane. Esterification of carboxylic acids and alcohols is a typical example of an equilibrium-limited reaction that produces ester and water as products. Limited by the equilibrium conversion, the aforementioned process generally shows low conversion, and an excess of reactants is required to enhance its value, determining an increase in the costs. Accordingly, removing reaction products from the reaction side, for instance, an ester or water, by pervaporation allows that, for the Le Chatelier principle, the reaction may be shifted toward the products. The esterification of LA with Et to produce EtLA and water is reported in the following [12][13][14][15][16][17][18][19][20]: Numerical models could be useful to avoid high experimental costs and to develop a better understanding of the effects of various parameters for the design and the fruitful application of PVMR in the levulinic acid esterification (LA-ESR) as well as for specific features and constraints such as the need of obtaining high conversion. To this purpose, the computational fluid dynamic (CFD) tool is a feasible method to simulate detailed liquid flow characteristics at any point of a membrane system. Indeed, the CFD approach can be used for the virtual prototyping of chemical reactors and separators. Since it is based on controlled volume methodology, the local variations of the fluid, thermal, and mass transport properties can be visualized in comparison to simple models and used to design a PVMR [21][22][23][24]. Nevertheless, to our best knowledge, there are not CFD-based studies about the application of PVMR for the LA-ESR. In this work, the LA-ESR reaction in a multi-bed PVMR (MB-PVMR) was theoretically analyzed, using the CFD method in order to evaluate the effects of the most significant operating parameters such as the reaction temperature, catalyst loading, feed molar ratio, and feed flow rate on the performance of the MB-PVMR in terms of LA conversion and water removal. Furthermore, the MB-PVMR performance was compared with those of a multi-bed traditional reactor (MB-TR) and discussed.

Developing of CFD Model
A CFD-based model was developed to simulate the performance of an MB-PVMR housing a poly(vinyl alcohol) (PVA) membrane, considering a two-dimensional and axisymmetric geometry. Figure 1 displays a simple scheme of the MB-PVMR (including four catalytic beds) adopted for EtLA production during LA-ESR reaction over a Smopex-101 catalyst. The main assumptions performed in the CFD simulations are summarized below: The film transport resistance at the interface of feed/membrane was considered negligible; • Absence of concentration polarization effects; • Membrane and catalyst not affected by deactivation.
Membranes 2021, 11, x FOR PEER REVIEW 3 of 14 the effects of the most significant operating parameters such as the reaction temperature, catalyst loading, feed molar ratio, and feed flow rate on the performance of the MB-PVMR in terms of LA conversion and water removal. Furthermore, the MB-PVMR performance was compared with those of a multi-bed traditional reactor (MB-TR) and discussed.

Developing of CFD Model
A CFD-based model was developed to simulate the performance of an MB-PVMR housing a poly(vinyl alcohol) (PVA) membrane, considering a two-dimensional and axisymmetric geometry. Figure 1 displays a simple scheme of the MB-PVMR (including four catalytic beds) adopted for EtLA production during LA-ESR reaction over a Smopex-101 catalyst. The main assumptions performed in the CFD simulations are summarized below: • Isothermal condition; • Steady-state regime; • The film transport resistance at the interface of feed/membrane was considered negligible; • Absence of concentration polarization effects; • Membrane and catalyst not affected by deactivation.

Governing Physical Equations
The governing equations in nonporous and porous zones are given in Table 1. The source term Si in the porous momentum equation is expressed as [25]: where μ is the liquid mixture viscosity, Kbr is permeability, and β F is the Forchheimer coefficient for packed-bed particles. Kbr and β F are defined as: The effective diffusion in porous media, D ei , depends on the structure of the porous material and the phases involved. In saturated or partially saturated porous media, the effective diffusivity is defined as: The fluid tortuosity for the Millington and Quirk model is:

Governing Physical Equations
The governing equations in nonporous and porous zones are given in Table 1. The source term S i in the porous momentum equation is expressed as [25]: where µ is the liquid mixture viscosity, K br is permeability, and β F is the Forchheimer coefficient for packed-bed particles. K br and β F are defined as:  The effective diffusion in porous media, D ei , depends on the structure of the porous material and the phases involved. In saturated or partially saturated porous media, the effective diffusivity is defined as: The fluid tortuosity for the Millington and Quirk model is: where ε, ρ, d cat , and D Li are the porosity, density of fluid, diameter of catalyst particle, and single-phase diffusion coefficients for the species in fluid, respectively. The porosity and particle diameter were set at 0.3124 and 0.01 mm, respectively. Porosity for the nonporous zone was set at 1. The flow pattern at the nonporous zone momentum equation is laminar flow. S i in the porous species equation is the sink/source term of component i, which accounts for the addition or removal of component i into the system for permeation through the membrane. In this work, since only water permeates from the retentate to the permeate side, it appears as a sink term on the retentate side and a source term on the permeate side. In other words, S i = 0 is considered for all components except for H 2 O, which is calculated as: where A is the membrane surface, V is the computational cell volume, M H 2 O is the water molar weight, and J H 2 O is the water permeating flux across the PVA membrane, calculated by Equation (8): The water permeating flux (J 0,H 2 O ) is reported at different temperatures in Table 2 [25].

Chemical Kinetic Reactions
The empirical reaction rate of LA-ESR reaction, conducted by Russo et al. [26,27] using Smopex-101 heterogeneous catalyst and the report self-catalyzed of this reaction is used as follows: r LA = r Et = −Rt and r LAEt = r H2O = Rt where R t and R 1 are the reaction rate at the catalytic and non-catalytic beds, respectively. Moreover, rate constants are defined according to Equation (12).

of 14
In Equation (10), ρ cat is the catalyst density and K the equilibrium constant expressed as: The parameters used of reaction rate constants are reported in Table 3. Table 3. Value of rate constant parameters used in the simulations of LA-ESR reaction [26,27].

Chemical-Physical Properties
The dependence of the liquid density on the temperature for each single component was evaluated according to the empirical expression reported below [27]: The coefficients of Equation (14) for each component are reported in Table 4.
where V i represents the molar volume at the normal boiling point (cm 3 /mol), calculated through the inverse of Equation (14), and µ mix is the viscosity of the reaction mixture (cP).
The term φ is defined as association factor and depends on the nature of the chemical component (φ LA = φ EtLA = 1.0, φ Et = 1.5 and φ H 2 O = 2.6). The φM i was therefore obtained for each component according to the Equation (16): where x j i is a molar fraction, in each step, determined for three components, using the i-th component and varying the j factor [28]. The temperature dependence on the viscosity was The coefficients of the Equation (17) for each component are reported in Table 5.

Boundary Conditions and Post Processing Definitions
In Table 6, the boundary conditions for retentate and permeate sides are shown. The following correlations were defined for describing the MB-PVMR performance in LA-ESR reaction: where LA in and LA out are the inlet and outlet levulinic acid molar flow rates, respectively, and H 2 O retentate and H 2 O permeate are the water molar flow rates in the retentate and permeate flow, respectively.

Numerical Method
Numerical simulations were performed using the commercial CFD package COMSOL Multiphysics 5.4, and the finite element method was used to solve the governing equations in the two-dimensional CFD model. Moreover, the pressure-velocity correction was performed using a coupling algorithm. Meanwhile, the equation solution was considered to be achieved when the residuals converged to values less than the magnitude of 10 −5 , and all the variable values were not changed with the iteration.

Mesh Independency
Another objective of the preliminary CFD simulations was to carry out a grid independence test. For this purpose, CFD simulations were carried out using different grid densities to find out the values beyond which the results become grid-independent. The investigated mesh numbers were 6020, 23,478, 36,240, 64,256, and 80,000. These simulations were carried out for MBPMR configuration during LA-ESR reaction (reaction temperature 343 K, feed flow rate 7 mm 3 /s, and molar ratio of Et/LA:1.25 and W = 8.6 g). The results for the grid independence tests are shown in Figure 2. The tracked parameter is the LA conversion. The results depict that beyond 48,000 meshes, the conversion shows a slowly decreasing trend, whereas it becomes constant at higher grid numbers. The finest grid was, therefore, identified as the grid density at which the solution becomes grid-independent. This grid number (48,000 meshes) was considered in all the final simulations discussed in the subsequent sections.  Figure 2. The tracked parameter is the LA conversion. The results depict that beyond 48,000 meshes, the conversion shows a slowly decreasing trend, whereas it becomes constant at higher grid numbers. The finest grid was, therefore, identified as the grid density at which the solution becomes grid-independent. This grid number (48,000 meshes) was considered in all the final simulations discussed in the subsequent sections.

Model Validation
The exactness of the CFD model was confirmed using trial information of Russo et al. [27]. Figure 3 illustrates the LA conversion versus the reactor length for the MB-TR during the levulinic acid esterification reaction. The working conditions for the TR model are T = 333 K, Et/LA = 1:1 and volume flow rate = 16.7 mm 3 /s, W = 8.6 g. As shown in the latter figure, the simulated and the experimental values are in suitable agreement.

Model Validation
The exactness of the CFD model was confirmed using trial information of Russo et al. [27]. Figure 3 illustrates the LA conversion versus the reactor length for the MB-TR during the levulinic acid esterification reaction. The working conditions for the TR model are T = 333 K, Et/LA = 1:1 and volume flow rate = 16.7 mm 3 /s, W = 8.6 g. As shown in the latter figure, the simulated and the experimental values are in suitable agreement.

Concentration and Velocity Distributions
The velocity distributions in the MB-PVMR, including the retentate and permeate zones, are shown in Figure 4b. The velocity values are indicated by the different colors of contours in agreement with the scale shown on the right side of the chart, where velocity distribution at

Concentration and Velocity Distributions
The velocity distributions in the MB-PVMR, including the retentate and permeate zones, are shown in Figure 4b. The velocity values are indicated by the different colors of contours in agreement with the scale shown on the right side of the chart, where velocity distribution at the nonporous zone follows a parabolic form in both sides of the MB-PVMR. Indeed, the maximum values of velocity are observed in the central area of the permeate and retentate sides. On the other hand, in the porous zone (catalyst bed), the velocity distribution is uniform in the catalyst bed. Figure 4a illustrated the velocity profile in porous and nonporous zones. Mole fraction distribution of H 2 O in both sides of the MB-PVMR is the most important parameter for the optimization of the process.

Concentration and Velocity Distributions
The velocity distributions in the MB-PVMR, including the retentate and permeate zones, are shown in Figure 4b. The velocity values are indicated by the different colors of contours in agreement with the scale shown on the right side of the chart, where velocity distribution at the nonporous zone follows a parabolic form in both sides of the MB-PVMR. Indeed, the maximum values of velocity are observed in the central area of the permeate and retentate sides. On the other hand, in the porous zone (catalyst bed), the velocity distribution is uniform in the catalyst bed. Figure 4a illustrated the velocity profile in porous and nonporous zones. Mole fraction distribution of H2O in both sides of the MB-PVMR is the most important parameter for the optimization of the process. Figure 5 sketches the distribution of H2O and levulinic acid concentrations as well as the permeation of water by membrane during production. The mass transfer flux in the shell and tube sides is composed of diffusion and convection. The convection mechanism tends to transfer H2O to the outlet of MB-PVMR owing to the high contribution of velocity in the vertical direction. The H2O is transferred toward the membrane by diffusional mass transfer that is responsible for the H2O removal.

Evaluation of Operating Parameters Effects
As reported in Table 7, after the preliminary CFD analysis, simulations were carried out to understand the effects of the most important operating parameters on the performance of the dense MB-PVMR in terms of LA conversion and water removal during LA-ESR reaction. The operating parameters were firstly checked in the model validation and, successively, we simulated the membrane reactor analyzing the mutual influence of them each other. This approach made it possible to understand from a theoretical point of view which kind of combination of operating conditions may determine better performance in terms of simulated conversion and water removal. In particular, the simulation sets can be subdivided into four parts, in which reaction temperature, catalyst loading, feed molar ratio, and feed flow rate values were changed for four cases of MB-PVMR and MB-TR.

Evaluation of Operating Parameters Effects
As reported in Table 7, after the preliminary CFD analysis, simulations were carried out to understand the effects of the most important operating parameters on the performance of the dense MB-PVMR in terms of LA conversion and water removal during LA-ESR reaction. The operating parameters were firstly checked in the model validation and, successively, we simulated the membrane reactor analyzing the mutual influence of them each other. This approach made it possible to understand from a theoretical point of view which kind of combination of operating conditions may determine better performance in terms of simulated conversion and water removal. In particular, the simulation sets can be subdivided into four parts, in which reaction temperature, catalyst loading, feed molar ratio, and feed flow rate values were changed for four cases of MB-PVMR and MB-TR. A further parameter that strongly affects the MB-PVMR performance is the feed molar ratio (Et/LA). High conversion of LA in higher content of Et is reachable. As shown in Figure 6, by increasing the feed molar ratio, LA conversion is enhanced. According to the LA-ESR reaction, an increase in Et/LA ratio may shift the reaction toward the products. Therefore, according to Le Chatelier's principle, the LA conversion and water production can be improved at a higher Et/LA value. In the meantime, H 2 O removal is also improved ( Figure 6). By considering the H 2 O compositions in both permeate and retentate sides, the H 2 O permeation driving force is increased at a higher feed molar ratio.

Effect of Feed Flow Rate on PVMR Performance
Residence time can be directly responsible for the conversion variation. It is clear that a higher flow rate may be responsible for a reduction in the residence time. As illustrated in Figure 7, an increase in feed flow rate reduced LA conversion. Similar to the previous case, LA conversion in the MB-PVMR is higher than that of the MB-TR. The higher the amount of H2O in the reaction zone, the higher the membrane surface required to overcome the reduction in H2O removal. Therefore, under a constant membrane area, the H2O removal is depleted. So, a feed flow rate increase from 2 to 10 mm 3 /s reduced the LA conversion from ~ 92% to ~ 78% in the MB-PVMR and from ~ 72.8% to 72.26% in the MB-TR.

Effect of Feed Flow Rate on PVMR Performance
Residence time can be directly responsible for the conversion variation. It is clear that a higher flow rate may be responsible for a reduction in the residence time. As illustrated in Figure 7, an increase in feed flow rate reduced LA conversion. Similar to the previous case, LA conversion in the MB-PVMR is higher than that of the MB-TR. The higher the amount of H 2 O in the reaction zone, the higher the membrane surface required to overcome the reduction in H 2 O removal. Therefore, under a constant membrane area, the H 2 O removal is depleted. So, a feed flow rate increase from 2 to 10 mm 3 /s reduced the LA conversion from~92% to~78% in the MB-PVMR and from~72.8% to 72.26% in the MB-TR. a higher flow rate may be responsible for a reduction in the residence time. As illustrated in Figure 7, an increase in feed flow rate reduced LA conversion. Similar to the previous case, LA conversion in the MB-PVMR is higher than that of the MB-TR. The higher the amount of H2O in the reaction zone, the higher the membrane surface required to overcome the reduction in H2O removal. Therefore, under a constant membrane area, the H2O removal is depleted. So, a feed flow rate increase from 2 to 10 mm 3 /s reduced the LA conversion from ~ 92% to ~ 78% in the MB-PVMR and from ~ 72.8% to 72.26% in the MB-TR.

Effect of Catalyst Loading on PVMR Performance
The catalyst performance during the reaction process is restricted by the catalyst load. Figure 8 depicts the impact of the catalyst load on the LA conversion and H 2 O removal in the MB-TR and MB-PVMR. Looking at the profiles of Figure 8, it is possible to observe that the increase in the catalyst load promotes the LA conversion. A slightly higher LA conversion in the MB-PVMR than that in the MB-TR is observed. The reason for that is attributed to the water separation in the MB-PVMR. In fact, by increasing the catalyst load, the conversion is enhanced and, therefore, H 2 O production is improved, leading to an increase in the driving force for a larger H 2 O removal. The catalyst performance during the reaction process is restricted by the catalyst load. Figure 8 depicts the impact of the catalyst load on the LA conversion and H2O removal in the MB-TR and MB-PVMR. Looking at the profiles of Figure 8, it is possible to observe that the increase in the catalyst load promotes the LA conversion. A slightly higher LA conversion in the MB-PVMR than that in the MB-TR is observed. The reason for that is attributed to the water separation in the MB-PVMR. In fact, by increasing the catalyst load, the conversion is enhanced and, therefore, H2O production is improved, leading to an increase in the driving force for a larger H2O removal. On the basis of reaction rate (Equation (11)) and H2O permeation flux (Table 1), which depends on the temperature through an Arrhenius-type equation, by raising the temperature, both LA conversion and H2O permeation flux increase. The prediction of CFD model-

Effect of Temperature Reaction on PVMR Performance
On the basis of reaction rate (Equation (11)) and H 2 O permeation flux (Table 1), which depends on the temperature through an Arrhenius-type equation, by raising the temperature, both LA conversion and H 2 O permeation flux increase. The prediction of CFD modeling in both MB-TR and MB-PMR as a function of the reaction temperature is shown in Figure 9. LA conversion increases from 77.5% to 86% and H2O removal from 77.5% to 79% as the reaction temperature raises from 323 to 363 K. The increase in the H2O removal reveals that the effect of temperature on H2O permeation is more relevant than its effect on water production. On the other hand, due to a higher H2O removal at higher reaction temperatures, the difference between LA conversion in MB-PVMR and MB-TR results to be consequently larger.

Conclusions
In this theoretical work, LA-ESR reaction by PVMR use was studied for producing ethyl levulinate as a significant fuel additive. An isothermal and 2D CFD-based model was adopted to investigate the performance of multi-bed PVMR configuration in terms of LA conversion and H2O removal, comparing the results with an equivalent multi-bed TR. The simulation analysis was particularly focused on the effects of several parameters such as reaction temperature, feed molar ratio, feed flow rate, and catalyst loading on the reactors, and the simulations well matched the experimental literature data. The most important aspect of this work was that, under all the considered operating conditions, the performance of the MB-PVMR was always better than those of the equivalent MB-TR. Furthermore, the best simulation results were reached at 343 K, 2 bar, and feed molar ratio equal to 3, obtaining a levulinic acid conversion equal to 95.3% and 91.1% of water removal. As a further result, the CFD tool was capable of predicting the optimal conditions to operate the MB-PVMR in this work. It pointed out also the benefits and drawbacks of the MB-PVMR during ESR reaction, demonstrating how it may constitute a valid strategy for the scale-up of this technology and the future commercialization of PVMR systems in added-value ester production processes for developing fuel additive industries.  LA conversion increases from 77.5% to 86% and H 2 O removal from 77.5% to 79% as the reaction temperature raises from 323 to 363 K. The increase in the H 2 O removal reveals that the effect of temperature on H 2 O permeation is more relevant than its effect on water production. On the other hand, due to a higher H 2 O removal at higher reaction temperatures, the difference between LA conversion in MB-PVMR and MB-TR results to be consequently larger.

Conclusions
In this theoretical work, LA-ESR reaction by PVMR use was studied for producing ethyl levulinate as a significant fuel additive. An isothermal and 2D CFD-based model was adopted to investigate the performance of multi-bed PVMR configuration in terms of LA conversion and H 2 O removal, comparing the results with an equivalent multi-bed TR. The simulation analysis was particularly focused on the effects of several parameters such as reaction temperature, feed molar ratio, feed flow rate, and catalyst loading on the reactors, and the simulations well matched the experimental literature data. The most important aspect of this work was that, under all the considered operating conditions, the performance of the MB-PVMR was always better than those of the equivalent MB-TR. Furthermore, the best simulation results were reached at 343 K, 2 bar, and feed molar ratio equal to 3, obtaining a levulinic acid conversion equal to 95.3% and 91.1% of water removal. As a further result, the CFD tool was capable of predicting the optimal conditions to operate the MB-PVMR in this work. It pointed out also the benefits and drawbacks of the MB-PVMR during ESR reaction, demonstrating how it may constitute a valid strategy for the scale-up of this technology and the future commercialization of PVMR systems in added-value ester production processes for developing fuel additive industries.

Conflicts of Interest:
The authors declare no conflict of interest.