Polyhydroxybutyrate Production from Natural Gas in A Bubble Column Bioreactor: Simulation Using COMSOL

In this study, the simulation of microorganism ability for the production of poly-β-hydroxybutyrate (PHB) from natural gas (as a carbon source) was carried out. Based on the Taguchi algorithm, the optimum situations for PHB production from natural gas in the columnar bubble reactor with 30 cm length and 1.5 cm diameter at a temperature of 32 °C was evaluated. So, the volume ratio of air to methane of 50:50 was calculated. The simulation was carried out by COMSOL software with two-dimensional symmetric mode. Mass transfer, momentum, density-time, and density-place were investigated. The maximum production of biomass concentration reached was 1.63 g/L, which shows a 10% difference in contrast to the number of experimental results. Furthermore, the consequence of inlet gas rate on concentration and gas hold up was investigated Andres the simulation results were confirmed to experimental results with less than 20% error.


Introduction
The facility of use and desirability of plastic materials properties led to their growing utilization in packaging and food industries [1]. At present, 30% of urban waste includes plastic waste. The extremely long durability of plastic (about three hundred years) has led to many environmental problems and destroyed the attractiveness of cities and nature. Burning polymeric lesions causes air pollution and according to the type of material used in their preparation, gases, such as hydrogen cyanide and other hazardous gases are released into the environment [1,2]. On the other hand, recycling has a lot of economic problems and costs. Since 1970, with the decline of the landfill problem worldwide, the issue of the use of biodegradable polymers was raised [3]. Polymers can be generally classified into two major biodegradable and non-degradable groups [4]. Non-degradable polymers, such as polyethylene, polypropylene, and polystyrene are produced from monomeric sources of oil and are of PHB is simulated, which precipitates the production process and reduces costs using the results of the experiment.

Materials and Methods
In order to examine the simulation of the process by using COMSOL software (Version 5.2: downloady.ir) and the production review in the bioreactor, first, a number of effective simplifying assumptions were selected: • The temperature of the system is always constant at 32 • C.

•
The physical properties of the solute with time are discarded.

•
The velocity of all compounds in the same phase is equal.

•
The culture medium and the microorganism mixture are considered as a single phase.

Reaction Kinetics
The production of PHB involves a type of methanotrophic bacterium, such as M. hirsuta, which is an aerobic and gram-negative bacteria that can produce PHB from the serum pathway. This methylotrophic bacterium of type II is used to investigate the production of PHB in the bioreactor, and the mathematical model equations provided by Equation (1) [35]: while C PHB is the dry weight of cell (g/L), K d is the cell death rate (s/1), and µ is specific growth rates for cell mass (1/s) (Appendix A). The modular kinetic model for producing PHB by Equation (2): In Equation (2), µ max is the maximum specific growth rate of microorganisms (1/s), S is the concentration of limited substrate for growth (g/L), and K S is the Monod constant (g/L). The kinetic of the Monod causes the substrate outlet concentration to be low in the CSTR because the K S is usually small. In order to create a delay phase, it is suggested that K d changes with time in Equation (3): α is the reversal of fixed cell death time (s −1 ), K d (∞) is the infinite cell death rate (s −1 ), and t is the time in seconds. It should be noted that for the implementation of the simulation reaction, the following simple reaction to the metabolic reactions is used for simulation:

The Equation of the Governing Model
The bubble flow model with two phases (diffused air in the liquid phase) was used in this study. In this model, the equations of continuity, momentum, and energy are solved for each step. The equation of motion is calculated by Equation (4) [34,35]: In this equation, u l indicates the velocity value in (m/s), p is the pressure in (Pa), and φ is the volume fraction indicated with m 3 /m 3 . ρ is density value with kg/m 3 , g is the gravity unit with m/s 2 , F is (N/m 3 ), µl is the dynamic velocity of the liquid replaced Pa.s in the equation. The values I and g, respectively, show the values of the liquid phase and the gas phase. The right-hand side of Equation (4) shows all forces that involve gradient, pressure, stress, adhesion, gravity, and force between the two phases, such as pulling, lifting, and virtual collective forces. In this study, a tension force is involved in the model. Based on the explanation, the continuity equation is written as Equation (5): (5) and the gas phase transfer equation is calculated as in Equation (6): m gl is the mass transfer rate from gas to liquid (kg/m 3 ). The gas velocity ug is equal to the sum of the velocity of Equation (7): U slip is the relative velocity among phases and U drift is the drift velocity. The physics relation calculates the density of gas from the ideal gas law by Equation (8): M is the molecular weight of the gas (kg/mol), R is the ideal gas constant (J/(mol·K) 3/3141472), and T is the temperature (K). p ref is a scalar variable being 1 at (1 at or 101.325 Pa) as a default. While a drift velocity is calculated by Equation (9) U drif = µ∇φ g ρφ g (9) µ is the effective viscosity, which causes it to fall. By putting Equations (9) and (7) in (6), we will have Equation (10): The drift velocity of the transfer equation is introduced in the gas transfer equation. This means that the gas transport equation is actually implemented in the physics interface. The equation of the bubble flow equation is relatively simple but it can indicate non-physical behavior. An artificial accumulation of bubbles, for example, is at the base of the walls in which the pressure gradients raise the bubbles while the bubbles have no place to go and there is no model for modifying the amount of gas fraction to grow. In order to prevent this, µ is set to µl for laminar. The only clear effect in most cases when the bubble flow equations are applicable is that the non-physical accumulation of bubbles reduces. The small effective viscosity in the transfer equation for ϕ g has beneficial effects on the numerical properties of the equation system.

Simulation Operations
In order to define the system, the properties and repercussions contained in COMSOL software were given to the simulation system. COMSOL simulation software has a complete set of information, properties, and constancy of materials, and if there is no specific information for a system, it can be found with the help Perry's handbook and enter into the physics of the problem. The geometry of the system was cylindrical with a radius of 0.015 m and a stanchion at the end of the cylinder with a radius of 0.1 mm. In order to prevent computing and convergence, and in view of the symmetry in the given geometry, the problem was solved in the symmetric two-dimensional mode, which did not have an effect on the overall solution due to the application of boundary conditions. All three components of water, methane, and ammonium were regarded as a phase. In this way, the liquid environment indicates the environment of the wastewater in which the ammonium and methane are soluble, and the properties of the bacteria are applied to the environment. The reaction occurred to produce biomass of an irreversible first-degree reaction with a reaction constant k 1 = 10 −5 in form of Va → W with the reaction rate, r = k 1 C va . The system temperature was set to 305.15 K. Since the environment is water, the incompressible flow with density and dynamic viscosity, respectively, were regarded as 1000 kg/m 3 and 10 -3 Pa.s. The initial pressure level (P 0 ) was considered zero and the permeation coefficient of the environment 10 -9 m 2 /s applied. The system under the boundary conditions of sleep with the equation u = 0 means that the system boundary was fixed and the simulation conducted in the symmetric two-dimensional mode. Figure 1 shows a schematic view of system geometry in the two-dimensional mode.

Resolution Independence Analysis of Numerical Grid
The number of model components is designed to meet the conditions of numerical resolution independence from the computational grid size and gradually increases the number of components. The resolution independence test of the computational grid to a point in which the difference in response is less than 5%. The number of components in resolution independence tests of the numerical grid is given in Table 1.

Meshing
At this stage, the geometry of the system was considered under meshing. In a two-dimensional model, the symmetrical mesh arrangement was selected based on the free triangular model and fine mesh size. The geometry of this system was split into 34,279 components. In Figure 2, a schematic view of the system meshing is presented in a two-dimensional mode.

Resolution Independence Analysis of Numerical Grid
The number of model components is designed to meet the conditions of numerical resolution independence from the computational grid size and gradually increases the number of components. The resolution independence test of the computational grid to a point in which the difference in response is less than 5%. The number of components in resolution independence tests of the numerical grid is given in Table 1.

Meshing
At this stage, the geometry of the system was considered under meshing. In a two-dimensional model, the symmetrical mesh arrangement was selected based on the free triangular model and fine mesh size. The geometry of this system was split into 34,279 components. In Figure 2, a schematic view of the system meshing is presented in a two-dimensional mode.

Meshing
At this stage, the geometry of the system was considered under meshing. In a two-dimensional model, the symmetrical mesh arrangement was selected based on the free triangular model and fine mesh size. The geometry of this system was split into 34,279 components. In Figure 2, a schematic view of the system meshing is presented in a two-dimensional mode.

Results
In the production of PHB from natural gas, the objective of optimizing the composition of the culture medium is to provide important and sensitive food, increase the yield of the product, prevent product degradation, and decrease the formation of harmful side effects.

The Results of Resolution Independence Analysis of Numerical Grid
As indicated in Table 2, three large, medium, and fine grid sizes were applied separately in the model and the final concentration of the biomass in each mesh size was studied to examine the sensitivity of the computational grid. Table 3 indicates the final concentration of biomass in the last step of solving equations in three computational grids.  Table 3. Final concentration of biomass in the last step of solving equations in three computational grids.

Concentration Contour
Biomass (C 5 H 7 NO 2 ) is generated under the applied conditions according to the reaction from Va → W. After the simulation, the concentration contour was determined at 14 min. Figure 3 displays the concentration contour of biomass. Figure 3 shows the concentration gradient is from the center toward the wall. As indicated in Figure 4, the concentration contour at time t = 2 min and t =10 min can be seen for the circulation and mixing of PHB [36].
Biomass (C5H7NO2) is generated under the applied conditions according to the reaction from Va → W. After the simulation, the concentration contour was determined at 14 min. Figure 3 displays the concentration contour of biomass. Figure 3 shows the concentration gradient is from the center toward the wall. As indicated in Figure 4, the concentration contour at time t = 2 min and t =10 min can be seen for the circulation and mixing of PHB [36].

Concentration Variations versus Time
The changes in the concentration of biomass were studied with time and the results are shown in Figure 5. As can be seen in this curve, the final concentration of biomass is almost 14 mol/m 3 , which is less than 5% in comparison to laboratory values (14.5-mol/m 3 ). Figure 6 indicates concentrationlocation changes in the two-dimensional model. As can be observed, the concentration at the beginning of the reactor is maximized because of the rapid reaction.

Concentration Variations versus Time
The changes in the concentration of biomass were studied with time and the results are shown in Figure 5. As can be seen in this curve, the final concentration of biomass is almost 14 mol/m 3 , which is less than 5% in comparison to laboratory values (14.5-mol/m 3 ). Figure 6 indicates concentrationlocation changes in the two-dimensional model. As can be observed, the concentration at the beginning of the reactor is maximized because of the rapid reaction.

Concentration Variations versus Time
The changes in the concentration of biomass were studied with time and the results are shown in Figure 5. As can be seen in this curve, the final concentration of biomass is almost 14 mol/m 3 , which is less than 5% in comparison to laboratory values (14.5-mol/m 3 ). Figure 6 indicates concentration-location changes in the two-dimensional model. As can be observed, the concentration at the beginning of the reactor is maximized because of the rapid reaction.

Concentration Variations versus Time
The changes in the concentration of biomass were studied with time and the results are shown in Figure 5. As can be seen in this curve, the final concentration of biomass is almost 14 mol/m 3 , which is less than 5% in comparison to laboratory values (14.5-mol/m 3 ). Figure 6 indicates concentrationlocation changes in the two-dimensional model. As can be observed, the concentration at the beginning of the reactor is maximized because of the rapid reaction.   Figure 7a,b indicate the results of spatial changes of fluid velocity in the bioreactor. On the right side of Figure 7, the velocity assigned to each color is represented.

Analysis of Variations in the Input Gas Velocity
The gas flow rate inside the reactor will influence the number of bubbles and the mass transfer  Figure 7a,b indicate the results of spatial changes of fluid velocity in the bioreactor. On the right side of Figure 7, the velocity assigned to each color is represented.   Figure 7a,b indicate the results of spatial changes of fluid velocity in the bioreactor. On the right side of Figure 7, the velocity assigned to each color is represented.

Analysis of Variations in the Input Gas Velocity
The gas flow rate inside the reactor will influence the number of bubbles and the mass transfer rate. Thus, the simulation was studied at four different velocities of 0.0015, 0.065, and 0.15 m/s. The results are displayed in Figure 8a-c.

Analysis of Variations in the Input Gas Velocity
The gas flow rate inside the reactor will influence the number of bubbles and the mass transfer rate. Thus, the simulation was studied at four different velocities of 0.0015, 0.065, and 0.15 m/s. The results are displayed in Figure 8a

Effect of Changing the Bubble Diameter on the Concentration
In Figure 9a,b, the change of concentration of the biomass can be observed with a bubble diameter of 3.5 and 1.5 mm.

Effect of Changing the Bubble Diameter on the Concentration
In Figure 9a,b, the change of concentration of the biomass can be observed with a bubble diameter of 3.5 and 1.5 mm.

Effect of Changing the Bubble Diameter on the Concentration
In Figure 9a,b, the change of concentration of the biomass can be observed with a bubble diameter of 3.5 and 1.5 mm.

Pressure Analysis
The relation between the amount of gravity and the bubbles are directed which are effected on pressure changes as presented in Figure 10.

Pressure Analysis
The relation between the amount of gravity and the bubbles are directed which are effected on pressure changes as presented in Figure 10.

Gas Accumulation
Gas accumulation is a dimensionless key dimensional parameter for design purposes identifying the phenomenon of transition in bubble column systems (Equation (11)). It is basically defined as the volume fraction of gas-phase occupied by gas bubbles. Similarly, liquid and solid phases can be determined as a liquid and solid phase coefficient. All studies examined gas accumulation because it plays an essential role in the design and analysis of bubble columns. The volume of gas accumulation with the mathematical relation was studied. In similar, the results achieved in the simulation show the very low gas accumulation in the lower part of the reactor and around the central axis, which can be observed in Figure 11a,b. Along the reaction, gas accumulation increases throughout the reactor but it is false since the accumulated gas is eliminated from the upper end of the reactor and ultimately the amount of accumulated gas at the bottom of the reactor is computed. The graph of gas volume coefficient throughout the reactor can be seen in Figure 11c, such as the value of gas volume coefficient and their increase and decrease tendency with time and position.

Gas Accumulation
Gas accumulation is a dimensionless key dimensional parameter for design purposes identifying the phenomenon of transition in bubble column systems (Equation (11)). It is basically defined as the volume fraction of gas-phase occupied by gas bubbles. Similarly, liquid and solid phases can be determined as a liquid and solid phase coefficient. All studies examined gas accumulation because it plays an essential role in the design and analysis of bubble columns. The volume of gas accumulation with the mathematical relation was studied. In similar, the results achieved in the simulation show the very low gas accumulation in the lower part of the reactor and around the central axis, which can be observed in Figure 11a,b. Along the reaction, gas accumulation increases throughout the reactor but it is false since the accumulated gas is eliminated from the upper end of the reactor and ultimately the amount of accumulated gas at the bottom of the reactor is computed. The graph of gas volume coefficient throughout the reactor can be seen in Figure 11c, such as the value of gas volume coefficient and their increase and decrease tendency with time and position. ε_g=u_g/(0.3+2u_g) (11) around the central axis, which can be observed in Figure 11a,b. Along the reaction, gas accumulation increases throughout the reactor but it is false since the accumulated gas is eliminated from the upper end of the reactor and ultimately the amount of accumulated gas at the bottom of the reactor is computed. The graph of gas volume coefficient throughout the reactor can be seen in Figure 11c

Shear Stress
In the present project, the shear stress in different input flows was studied indicating (20) and (21) the stress rate, which is based on the scale available on the right of these shapes, the results achieved in line with the expectation of the sign, which contributes to an increase in the stress rate at the vicinity of the oven in line with the increase in the rate of intake gas and discharge. Figure 12 indicates the stress rate over time.

Shear Stress
In the present project, the shear stress in different input flows was studied indicating (20) and (21) the stress rate, which is based on the scale available on the right of these shapes, the results achieved in line with the expectation of the sign, which contributes to an increase in the stress rate at the vicinity of the oven in line with the increase in the rate of intake gas and discharge. Figure 12 indicates the stress rate over time.
In the present project, the shear stress in different input flows was studied indicating (20) and (21) the stress rate, which is based on the scale available on the right of these shapes, the results achieved in line with the expectation of the sign, which contributes to an increase in the stress rate at the vicinity of the oven in line with the increase in the rate of intake gas and discharge. Figure 12 indicates the stress rate over time.

The Results of Resolution Independence Analysis of Numerical Grid
As Table 2 shows, the variance in the three large, medium, and fine grid sizes is less than 1% indicating the independence of the obtained responses on the type and size of the grid. Figure 3 shows that from the center toward the wall, concentrations are reducing being observed in Figure 3b. Thus, the largest value of production is located near the central axis. By enhancing the velocity of the areas around the center of the bioreactor, the maximum concentration, and with a distance less than the center of the percentage of production means that the increase in velocity would lead to a shorter spatial dimension to the final production.

Concentration Contour
Mousavi et al. examined the simulation of PHB production in an aircraft bioreactor with Fluent software. The average molecular concentration of PHB in the bioreactor was almost 0.331 and 0.3337 mol/L [36]. Figure 5 shows the changes in the concentration of biomass were investigated vs. time. The final concentration of biomass is 14 mol/m 3 . The researchers showed the production of PHB was 8 g/L in a bubble column bioreactor at optimized conditions [17]. In addition, Shah Hosseini et al. applied a dynamic optimization program written in MATLAB software to specify the optimal amount of carbon and nitrogen sources on PHB production. The proposed model presented was confirmed on the experimental data. Using the preferred feed strategy, PHB was enhanced by 100% [37]. The Euler-Euler approach is used to clarify two-phase motion equations that correspond slightly to the experimental data for both mid-velocity and motor velocities. LES provides better matching with simulation data through the k-ε model [38].

Velocity Contour
Mousavi et al. studied the fluid velocity fluctuations of the aircraft bioreactor in which the maximum fluid velocity can be seen in the input aperture at the bottom of the reactor. Similarly, the velocity of the fluid in the Down Comer is much lower than the Riser section [36].
Liquid velocity can directly effect production and methane gas enhances the growth structure of the biomass because these two parameters change in parallel. The flow rate can also change the behavior of the bioreactor and can change the production speed according to the size of the meshes designed [39,40].

Analysis of Variations in the Input Gas Velocity
The result of change of the gas flow rate inside the reactor vs. the number of bubbles and mass transfer rate (at four different velocities of 0.055, 0.015, 0.065, and 0.15 m/s) are displayed in Figure 8. As can be observed, as the gas flow rate increases, the content of biomass increases because of the increase in the number of bubbles and the number of encounters, which finally leads to an increase in the mass transfer rate and biomass production. The effects of gas and liquid velocity in a bubble pillar biomass with COMSOL software varying in the gas-liquid bubble column system, the velocity of gas and liquid with time and space in the column. Weather velocity vectors were achieved after a semi-stable mode at an input velocity of 0.001 m/s [41]. Such vectors indicated the velocity of the path along the path, which is useful in specifying the patterns of flow in the bubble column system. That contour indicated a gas concentration of 0.001 m/s. Figure 9 shows that the concentration of the biomass is reduced by increasing the diameter of the bubble to reduce the flow rate. In addition, it enhances the final concentration of biomass by decreasing the diameter of the bubble.

Pressure Analysis
As it can be observed in Figure 10, the pressure across the reactor from bottom to top reduces because of the presence of gravity on the bubbles.

Gas Accumulation
The graph of changes in gas volume throughout the reactor can be seen in Figure 11. In 2010, Mousavi et al. studied gas cumulating in aircraft bioreactors at numerous aeration rates [36]. The results obtained from the simulation were compared to the results of different mathematical relations related to the bubble columns and aircraft reactors. The simulation results are relatively suitable in comparison to the relationships and it approves the reliability of the simulation. For the aeration rate of 30 L/g. min, the mean diameter of the bubble was 5.1 mm [36].
Gas accumulation occurs during the five different phases of aeration. As expected, the amount of gas accumulation increased with aeration. The highest volume is always increased initially and then gas accumulation decreases because of the gas outflow from the bioreactor. However, in cases with higher aeration rates (40 and 50 L/min), the final amount of gas accumulation is higher than the peak value. In higher aeration conditions (40 and 50 L/min), the amount of gas accumulation is higher than the initial rate of gas accumulation in comparison to the lower aeration rate due to the accumulation of gas in the liquid phase.

Shear Stress
Wall shear stress is a critical parameter for determining energy transfer and movement in a two-phase flow system. The shear stress properties play an essential role in understanding the internal state of the two-phase flow because the liquid velocity and the highest grade of turbulence fluctuations are specified with the highest gradient. In order to clarify the mechanism of the two-phase bubble flow, it is significant to know the performance of the bubbles in the fluid flow. The movement of isolated bubbles seems to be related to the bubble distortion, the location of the injection point, and the average fluid flow rate. They found the only bubbles with diameters less than 5 mm along the wall while all the spherical bubbles and elliptical bubbles with a diameter greater than 5 milligram meters in the core of the stream [42,43]. Mousavi et al. in a study on shear stress in aircraft bioreactors specified the average shear stress in the bioreactor by applying the equations in FLUENT software for five different conditions. As expected, aeration enhances the shear stress in the bioreactor. Since airborne bioreactors and bubble columns fail at having a propellant, their tensions are lower than those of cohesive bioreactors, and the shear stress in the bioreactor is not likely to be significant, while the ability of the software to calculate shear stress is critical. The greatest shear stress relates to areas where gas is flowing from the opening. Vortices are located at the vicinity of upstream and downstream sections as well as the walls with higher shear stresses compared to other points [36].

Conclusions
In this study, a simulation of microorganisms capability for the production of PHB from natural gas (source of carbon) was carried out. The error interval is normally up to 20% for the difference between laboratory data and simulation because the resulting error can be because of simple simulation assumptions; thus, the simulation results were in good accordance with empirical results. Regarding the concentration contour, concentration is the lowest in the center of maximum concentration near the wall. The concentration gradient from the center is toward the wall and produces the highest percentage of production near the central axis of the bioreactor. Substrates and microorganisms with the environment were considered homogeneous, thus we would not have much difference in production compared to the place. Enhancing the gas flow rate has a direct relation to the production rate. Based on the three-velocity concentration-time graph, enhancing the intake rate of the agitator stirring increases and the percent of production enhances. The flow pattern in a bioreactor is dependent on the diameter of the bioreactor and the flow rate of the inlet gas. In this study, the bubble flow studied it placed due to the low velocity of the inlet gas (0.015 m/s) and the low diameter of the bioreactor (1.5 cm). Because of the concentration contour with enhancing the flux, the number of bubbles increases and results in more mass transfer and an increase in the value of biomass production. With increasing fluid velocity, a decrease can be seen in the amount of gas accumulation because liquid bubbles are rapidly displaced by liquid at higher velocities. Moreover, increasing fluid velocity decreases the duration of bubble stays.
Furthermore, based on the results in volume concentration, 50/50 from air and methane, the highest rate of microorganism growth and PHB production was achieved. According to presented graphs, such as time concentration and location concentration, the final amount of biomass production extracted was 1.6338 g/L under optimal conditions. The reason for selecting the volume of 50/50 was the coordination with the Taguchi algorithm that is based on research in the same conditions is the ideal conditions for evaluation of chosen sizes of the system. In addition, the amount of mesh size in the system was not significantly affected by the final concentration of production, so the ideal condition was regularly 1.6338 g/L. Thus, in different conditions of meshes, less than 1% difference was observed. In addition, the amount of gas accumulation decreases with increasing gas velocity. Acknowledgments: Authors wish to acknowledge all administrative and technical support supports given by Tehran University.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
List of symbols in equations: