CFD-PBE Modelling and Simulation of Enzymatic Hydrolysis of Cellulose in a Stirred Tank

Corresponding Author: Norazaliza Mohd Jamil Faculty of Industrial Sciences and Technology, Universiti Malaysia Pahang, Lebuhraya Tun Razak, 26300 Gambang, Kuantan, Pahang, Malaysia Email: norazaliza@ump.edu.my Abstract: Renewable energy or biofuel from lignocellulosic biomass is an alternative way to replace the depleting fossil fuels. The production cost can be reduced by increasing the concentration of biomass particles. However, lignocellulosic biomass is a suspension of natural fibres and processing at high solid concentration is a challenging task because it will affect the mixing quality between the enzyme and cellulose particles and the generation of sugars. Thus, understanding the factors that affect the rheology of biomass suspension is crucial in order to maximize the production at a minimum cost. Our aim was to develop a solution strategy for the modelling and simulation of biomass suspension during enzymatic hydrolysis. The complete model was solved using the DAE-QMOM technique in a finite-element software package, COMSOL. Essentially, we made a clear connection between the microscopic and macroscopic properties of biomass slurries undergoing enzymatic hydrolysis. The results showed that the quality of mixing within a reactor is crucial in optimizing the hydrolysis product. The model improved the predictive capabilities, hence increasing our understanding on the behaviour of biomass suspension.


Introduction
Enzymatic hydrolysis of cellulose is a complex phenomenon. It consists of several effects that occur simultaneously in a biomass suspension reactor, such as fragmentation of cellulose chains, in homogeneity in the tank, mixing effect and adsorption of enzymes to the cellulose substrates. In fact, even though numerous issues have been explored, they cannot be considered as fully covered.
It was discovered that advection did represent a significant phenomenon that could increase the number of cellulose particles generated during the hydrolysis process. Besides advection and the diffusion factor, another critical aspect in the enzymatic hydrolysis of cellulose for bioethanol production is the contact between the enzymes and the cellulose substrates (Jager and Buchs, 2012). The problems that arise at present are very low solubility of cellulose substrates and as for the economic issue; costly high enzyme concentrations are commonly applied.
Hence, this paper depicted the development of a solution strategy in resolving the modeling, as well as the simulation of the PBE-advection-diffusion system, coupled with a fluid medium surrounding the cellulose particles. The cellulose particles were suspended in a continuous fluid that was assumed to obey the incompressible Navier-Stokes equations. Besides, the mixing between biomass suspension and cellulase enzymes was simulated in a cylindrical stirred tank, mechanically agitated by a stirrer. In addition, the influence of the mixing speed of the stirred tank and different locations of enzyme injection points in the mixer to the hydrolysis yield were investigated numerically. Meanwhile, the model of Griggs et al. (2012a) for the hydrolysis of cellulosic biomass was extended and modified to accommodate the Computational Fluid Dynamics (CFD), which led to strong spatial in homogeneities in the flow and the species concentration fields. Furthermore, the capability of DAE-QMOM approach (Gimbun et al., 2009) was tested to solve the coupled model.
Many theoretical and experimental studies have been described in the literature, trying to capture the important aspects of the enzymatic hydrolysis process (Abnisa et al., 2011;Bansal et al., 2009;Carvajal et al., 2011;Fan and Lee, 1983;Griggs et al., 2012a;2012b;■■■ Himmel et al., 2007;Igarashi et al., 2009;Kleman-Leyer et al., 1996;Limayem and Ricke, 2012;Okazaki and Moo-Young, 1978;Selig et al., 2008;Srisodsuk et al., 1998). From these studies, several kinetic models based on discrete and continuous formulations on population balance equations had been proposed. However, to the best of our knowledge, the coupling of enzymatic hydrolysis kinetic model with fluid dynamic has not been reported. Nevertheless, the Standard Method of Moments (SMM) and Quadrature Method of Moments (QMOM) had been employed in conjunction with CFD to simulate the process of precipitation of barium sulfate (Cheng et al., 2009). In another study, Marchisio et al. (2003) implemented the QMOM in a CFD code FLUENT in modeling simultaneous aggregation and breakage problem.
This paper is organized as follows. First, the development of phase-space kinetic is illustrated for the model formulation. Next, in the third section, the setup of the 2-D stirred tank simulated in COMSOL is presented. Besides that, the mixing in the stirred tank is discussed and then, the results and discussion are provided by using the model developed in this study. Finally, a conclusion of this paper is given. This study is expected to provide a valuable insight concerning the features of the breakage phenomena during the enzymatic hydrolysis of cellulose in a stirred tank.

Phase-Space Kinetic Model
The population balance for the enzymatic hydrolysis of cellulose problem was originally expressed in terms of the number of the density function, p(χ, t) with the particle length χ as its internal coordinate (Griggs et al., 2012a). The PBE describes how the particle size distribution changes as time progresses due to polymer breakage during cellulosic hydrolysis. Furthermore, it is assumed that the fluid environment plays an explicit role in particle behavior. In this case, the distribution of particles depended not only on internal coordinates, but also on its location in the physical space, referred to as external coordinates, x, here. The external coordinate was utilized to locate the physical positions of suspended cellulosic particles, which were governed by flow advection and particle diffusion (Ramkrishna, 2000). From these coordinates, the PBE for an enzymatic hydrolysis process can be governed by: where, u is the mean velocity vector, D p and D pB are the particle diffusion constants for the length χ of p and p B , respectively. Here, p denotes the population distribution of enzyme-accessible cellulose chains of length χ, p B is the population distribution of CBH 1 -threaded cellulose chains of length χ and R is the radius of the cellulose particles. The source terms f p , f pB and f R represent the reaction terms due to the breakage process by EG 1 and CBH 1 enzymes of the cellulose chains, which are given respectively by: The equations in (1) are the integro-partial differential equations, which describe how the dispersed phase evolves in time, space and size coordinates. The problem apparently became more complicated when external coordinates appeared in the models. Hence, we nondimensionalized the system so as to reduce the number of parameters involved and to make the equations compact.
The Method of Moments (MOM) was first introduced by (Hulburt and Katz, 1964) and it forms the basis for many promising solution strategies. The basic idea of MOM is to transform the original PBE into a problem in identifying the moments in size distribution. However, applying the definitions of moments to the reduced model equations results in a closure problem for the system. Therefore, a practical approach to solve the closure problem is the QMOM, introduced by (McGraw, 1997). In QMOM, the integral terms containing the size of the distribution function are approximated by numerical quadratures in the same fashion as in numerical integration. With this approach, the moments can be expressed as: and: where, n q is the order of the quadrature formulation, ξ i and L i are the particle lengths, while w i and c i are the quadrature weights for p and p B , respectively. About 60 Accordingly, p (0) is the total number of enzyme accessible cellulose particles and p (1) is the total number of monomeric glucans for enzyme accessible cellulose particles. Subsequently, by combining the algebraic equations from the QMOM definitions and the partial differential equations from the moment transport equations, the Differential Algebraic Equations (DAE) system was obtained. This method is known as the DAE-QMOM method, proposed by (Gimbun et al., 2009). Therefore, the full system can be expressed as: The developed mathematical model consisted of three parts: (i) The hydrodynamic core, a model to solve the incompressible Navier-Stokes equations; (ii) transport equation for advection and diffusion problems, (iii) PBE to describe the size distribution of cellulose particles. The schematic view of how these models are interconnected is given in Fig. 1. The coupled phasespace kinetic model with hydrodynamic is visualized as a single phase laminar flow. The behavior of the fluid is described by the incompressible Navier-Stokes equation as portrayed below: where, u is the velocity vector, p is pressure, ρ gives the fluid's density and µ denotes the dynamic viscosity. The velocity of the laminar flow calculated from the Navier-Stokes equations was used to solve the moment transport equations by providing the fluid convective velocity. Meanwhile, the reaction terms in the transport equations consisted of quadrature weights and quadrature abscissas variables, which can be solved with nonlinear algebraic equations. The DAEQMOM technique that originally solved the pure PBE problem, developed by (Gimbun et al., 2009), had been tested in solving the coupled fluid dynamics and population balance kinetic models for the enzymatic hydrolysis of biomass.

Setup of 2-D Stirred Tank in COMSOL
The reactor tank simulated in this study was a horizontal stirred mixer. The cylindrical mixing tank is represented by a 2D schematic cross-section of a tank, as illustrated in Fig. 2. The mixer tank was stirred by fourblade impellers and had four equally-spaced baffles attached to the wall to enhance mixing. The impellers rotated clockwise at a speed of 1-9 RPM and the fluid in the mixer consisted of biomass suspension. The cellulose molecules were long chains of mono-saccharides and as the fluid was non-Newtonian and viscous; Carbopol with a density at 208 kg/m 3 and viscosity at 4000Pa⋅s was simulated as a model fluid to mimick the properties of biomass suspension.
The geometry of the mixer was divided into two parts: The rotating domain that formulated the Navier-Stokes equations located in the small circle covering the impeller; and the fixed material coordinate system was applied on the outside ring covering the baffles. Both the inner and the outer domains were linked to each other by an assembly boundary, i.e., the line circle between the impeller and the wall of the mixer. Here, the line circle had a flux continuity boundary condition and it transferred momentum to the fluid in the outer domain. Besides, a no-slip boundary condition was imposed to the reactor's fixed walls. This was consistent with no inflow and outflow from the reactor, which was expected. At the impeller, the convection followed the flow, so the impeller boundaries were set to a convective ux. Furthermore, the CFD simulation in COMSOL 4.3b employed the following setting: Type, 2D; and Model: Rotating Machinery Laminar Flow (incompressible flow), Transport of Diluted Species and Domain ODEs and DAEs. The Rotating Domains feature was applied to specify the rotational frequency and the direction that represented the effect of an impeller on the fluid. Moreover, a single-phase model with PBE-advectiondiffusion, coupled with hydrodynamics, was developed and solved with the DAE-QMOM approach within COMSOL. It solved a momentum balance to fully describe the stirred fluid. It also solved a material balance to describe how the concentration distribution was affected by the stirring. Meanwhile, the seeding of the enzyme was done at a point in the mixer, while the mixing between the enzyme and the fluid was monitored.

Mixing in a Stirred Tank
The impellers in a stirred tank reactor were used to generate flow and mixing within the reactor. Mixing in reactors is a challenging problem that has attracted numerous studies, including experimental and fundamental theories. The interaction of chemical reactions and fluid mechanics can cause significant complexity to a system (Madras and McCoy, 2004). The efficiency of a reactor is determined by the efficiency of the mixing within the reactor, as good mixing increases the production process. Besides, increased fundamental understanding of hydrodynamics and mixing in such reactor leads to a reduction in power usage.
Moreover, several studies investigating the modeling of fluid dynamic in a stirred tank have been carried out. For example, (Guha et al., 2006;Micale et al., 2000) accounted for flow description within the reactor, but ignored the kinetic effects of the reactor performance. On the other hand, Madras and McCoy (2004) studied the mixing effect in a stirred tank incorporating the population balance equations but did not take into account the flow description within the reactor. Therefore, in order to study the effect of hydrodynamics on the mixing behavior, decoupling the flow and kinetics of the system are required.
In the case of the enzymatic hydrolysis of cellulose, the enzymes in a reactor need a sufficient contact with the cellulose substrates before the fragmentation of cellulose chains can take place. Hence, mixing within the reactor is required. More importantly, mixing is necessary to avoid areas with high concentrations of hydrolysis products (cellobiose and glucose), which would inhibit the enzymes (Jorgensen et al., 2007). This is because; a good enzyme distribution is necessary to facilitate high biomass conversion rates (Roche et al., 2009). Consequently, the production of hydrolysis yield can be affected by the mixing. Hence, a proper understanding of the mixing effects in the cellulosic hydrolysis process can lead to the reduction in cost and increase the profitability of operation.
Apart from that, a biomass suspension, together with a solution of enzymes, was stirred in a reactor. In the simulation carried out in this study, the enzyme was injected at a point inside the stirred tank reactor. Then, the agitation of the impeller caused the enzyme to spread to a different place in the domain. Thus, the quantification of the mixing that varied with time in the stirred tank had been desirable. This was done by integrating the difference between the actual concentration and the mean concentration throughout the domain using the following formulation (Comsol, 2007): where, c is the actual concentration and c mean is the mean concentration. Based on Equation 12, the MixingValue gave a relative value of how well the mixing inside the tank was and how far it was from the ideal mixing (Comsol, 2007). MixingValue is a key variable, which provides a way to evaluate mixing effectiveness and to compare the mixing systems. It provides information at which the tank attains a specified degree of uniformity. For the present study, the mixing of carbopol and enzymes in a stirred tank was modeled and simulated by using the COMSOL software. The fluid flow and the population balance equations were coupled for the cellulosic hydrolysis process through the DAE-QMOM approach, in which it bridged the gap between CFD and the phenomenological model. The goal was to predict the influence of mixing on the production from the enzymatic hydrolysis of cellulose. In order to do so, the effects of the mixing speed and enzyme feeding location on the hydrolysis yield in the system had been identified.
The conventional evaluation mixing was done through experiments in a laboratory. This approach is usually expensive, time-consuming and difficult. In this regard, a valid model simulation offers a better alternative, whereby one can examine various parameters of the mixing process within a shorter time and at less expense. The purpose of the present study was to evaluate the capability of numerical simulation, i.e., DAE-QMOM, to predict different features of the mixing process for cellulosic hydrolysis. COMSOL 4.3b was used to solve PBE that operated in conjunction with its flow calculations to predict the behavior of the system. The two-dimensional cross-sectional reactor that was constantly stirred by the movement of an impeller had been simulated. The problem was indeed challenging since the effects of fluid flow, transport of species and reaction kinetics had to be included in the model simulation.
Here, the biomass suspension was considered as a single-phase laminar flow. Figure 2 displays the velocity distribution of carbopol in a stirred mixer calculated with the Navier-Stokes equation with density ρ = 208 kg/m 3 and viscosity µ = 4000Pa⋅s, as well as a mixing speed of 1 rpm. The red arrows illustrate the flow pattern of the fluid, while its length indicates if the fluid moves slowly or fast. As depicted in the figure, the fluid near the walls moved slower than the fluid in the center of the tank. In addition, the highest velocity magnitude was located at the tip of the mixer blades. Figure 3 shows that the number-averaged chain length, (1) , decreased with time during the hydrolysis process in the coupled PBE-advection-diffusion with the hydrodynamics model. This result indicated that the long cellulose chains were fragmented into shorter ones by the actions of enzymes EG 1 and CBH 1 , which had been consistent with the previous findings in the earlier chapters. Furthermore, several tests on the mixer pertaining to the hydrolysis process were carried out. First, the effect of mixing speed on the hydrolysis yield was predicted. Then, the model was tested to show if the feeding location of enzymes brought any difference to the production of sugars.

Effect of Mixing Speed
On top of that, the effect of mixing speed to the cellulosic hydrolysis yield in the stirred tank was investigated. The simulation for various mixing speeds i.e., 1, 3, 5, 7 and 10 rpm, was conducted for two different enzymes feeding locations. The feeding locations were in the outer and the inner circle of the stirred tank. Meanwhile, Fig. 4 portrays the concentration of glucose at 0.158 mol/m 3 for enzyme feeding in the outer circle, while 0.365 mol/m 3 for the inner circle after 3 h at different mixing speeds. No significant difference was found between the glucose concentrations at different mixing speeds of a stirred tank for both feeding locations after 3 h of cellulosic hydrolysis. The glucose production was, therefore, insignificantly affected by the mixing speed within the tested range.
In addition, during hydrolysis, stirring or shaking in a mixer tank is required to enhance the contact between enzyme and cellulose particles and to increase the mobility of enzymes throughout the entire volume (Jorgensen et al., 2007;Mais et al., 2002). Product inhibition by glucose and cellobiose at increased concentrations at high solids loading plays a role for the decreasing conversion (Kristensen et al., 2009). Therefore, stirring in a tank helps avoid the areas with high concentration of hydrolysis products (cellobiose and glucose) that can inhibit the action of enzymes (Jorgensen et al., 2007). However, high-stirring speed is not recommended as it may result in enzyme denaturation and deactivation by the high shear forces (Cao and Tan, 2004;Tengborg et al., 2001).
The correlation between stirring and shaking speed and the generation of hydrolysis products has been studied previously, but in an experimental way using shake asks or stirred tank reactor (Jorgensen et al., 2007;Ingesson et al., 2001). Besides, prior studies have noted that the mixing speed of a mixer did not significantly affect the generation of hydrolysis products (Jorgensen et al., 2007;Ingesson et al., 2001). Hence, the simulation results obtained in this study matches those observed in earlier experimental studies which indicated that the coupled DAE-QMOM approach with hydrodynamic is capable of predicting the hydrolysis process. A major constraint in the enzymatic hydrolysis of biomass for ethanol production is the cost of cellulase enzymes and any strategy that can bring down the production cost of cellulases can significantly reduce the cost of bio-ethanol (Sukumaran et al., 2009). The ability of DAE-QMOM to simulate the hydrolysis system rather than an experimental test is beneficial in avoiding production cost due to expensive enzymes, machines and energy.

Effect of Enzyme Feeding Location
Moreover, the influence of the enzymes feed locations to the mixing quality and its correlation to the production of sugars during hydrolysis process had been looked into. To do these, the mixing value and the hydrolysis yield were measured when the enzyme was injected at four different points inside the stirred tank. The COMSOL program package was used for the model simulation that incorporated reaction kinetics and fluid flow properties. In cellulosic hydrolysis process, several phenomena occurred simultaneously. First, the particles began to break and due to the velocity of the fluid, they moved clockwise in the reactor. The diffusion and the advection velocity had been accounted as the main factors that affected particle distribution in the tank.
Besides, it was known that the flow field in the stirred tank had two main regions: The one around the impeller, where most breakages occur and the outer circle of the tank. Using visualization, a significant vortex of the enzyme inside the tank was observed. The system can be observed from the snapshots of the concentration field at different feed injection locations as shown in Fig. 5 and 6. Figure 5 visualizes the case where the feeding enzymes were located in the outer circle. Since the velocity had been rather slow at the seeding point, the initial transport was almost isotropic from the seeding point. It clearly showed that the enzyme was still undispersed and was concentrated in the outer circle of the reactor.
On the other hand, Fig. 6 shows the enzyme distribution in the mixer when the feeding location was in the inner circle. This feeding was in the faster velocity field in the center of the mixer and was spread in an azimuthal direction. It only took a few minutes more for the enzyme to be almost homogeneously distributed in the mixer. The enzyme was strongly dispersed due to the mixing effect that the recirculation zones had on the reactor. To sum up, higher localized concentration was observed when feeding was at the outer circle, as opposed to being almost uniformly spread throughout the domain when the feeding was in the inner circle.
Meanwhile, Fig. 7 illustrates a sketch depicting the differences in the feeding locations of enzymes in a stirred tank reactor, i.e., the (1) outer circle, (2) inner circle, (3) tip of the impeller and (4) on the impeller, whereas Fig. 8 shows the MixingValue as a function of time for different injection points of enzymes in the mixer. MixingValue is the difference of the actual concentration and the average concentration integrated over the domain. The increase at early times had been due to the inlet of the enzymes.  (3) and (4) yielded similar results. It was observed that the MixingValue curve for enzyme feed on the impeller had a peak above both the outer and the inner circle. Besides, the required time to reach equilibrium or homogenous distribution throughout the domain had been the longest for the enzyme feed on the impeller than in the outer or inner circle. In other words, the slowest spreading occurred when the feed injection was located on the impeller blades. Another essential point was that the enzyme feed injection location located in the inner circle had the lowest peak and the shortest time to reach equilibrium compared to the other two locations.
In a similar way, the evolution of glucose concentration corresponding to the feed injection location is shown in Fig. 9. The enzyme injection on the impeller yielded the least glucose concentration, followed by the outer circle. The highest glucose concentration was obtained by feeding the enzyme in the inner circle of the mixer. As time progressed, the enzyme spread more in the reactor and hydrolysis qualitatively increased in all cases. The glucose concentration increased over time as expected due to the reactions on the cellulose chains reacted by EG 1 and CBH 1 enzymes. Besides, this phenomena indicated a correlation between the MixingValue and the production of glucose for each location. In the enzymatic hydrolysis of the cellulose process, sugar productions go hand in hand with efficient mixing between the enzymes and the cellulose substrates. The outer circle was clearly worse than the inner circle feeding location in the mixing process, but better than the one on the impeller.   Lavenson et al. (2012) performed an experimental study using Magnetic Resonance Imaging (MRI), a cylindrical penetrometer and HPLC and reported that the spatial homogeneity in the distribution of enzymes (mixing quality) had a major influence on hydrolysis rate. Besides, according to (Hodge et al., 2008), mass transfer was a limiting factor in the decreased production of cellulase enzymes. Meanwhile, as explained by (Griggs et al., 2012b), predictive models for the enzymatic hydrolysis of cellulose enabled an improved understanding of the current limitations and might have pointed the way for a more efficient utilization of enzymes in order to achieve higher conversion yields.
On the other hand, the findings retrieved from this study revealed that the feeding location of enzyme in a stirred tank affected the mixing quality of the hydrolysis process. Therefore, the location of enzyme feeding must be chosen wisely to minimize diffusion limitations without denaturing the enzymes (at a high shear rate). Moreover, changing the feeding location from the outer ring and on the impeller to the inner circle increased the MixingValue. Hence, for a specific geometry of a stirred tank, the best location to feed the enzyme to the biomass suspension had been in the inner circle, close to the impeller, but not on the impeller.

Conclusion
In this study, the enzymatic hydrolysis of cellulose in a stirred tank reactor was modeled. The Griggs kinetic model was extended by adding diffusion, advection and hydrodynamics, which had a simultaneous effect on the system. The flow led was obtained by solving the singlephase Navier-Stokes equations with a standard laminar flow. The population balance equation was solved through the combination of Differential Algebraic Equations and Quadrature Method of Moments (DAE-QMOM). The DAE-QMOM was incorporated into an inhouse CFD code to simulate the enzymatic hydrolysis of cellulose in a stirred tank reactor. Besides, numerical simulation using 2-node QMOM indicated that the results showed a good qualitative agreement with the experimental results.
Hence, a DAE-QMOM approach was found to be effective for PBE-advection diffusion model with the incorporation of a CFD framework. This method had been shown to provide reasonable predictions at a reduced computation expense for a single-phase system. In addition, this approach enabled one to get a better understanding of the hydrolysis model as the owing fluid had been taken into account as the environment in the hydrolysis process. Besides that, further insight into the mixing process was obtained by predicting the effects of the mixing speed and enzyme feeding locations in the mixer. The quality of mixing within the reactor is crucial in optimizing the hydrolysis product. Moreover, the results from this paper showed a qualitative agreement with the experimental results obtained in previous studies. In conclusion, the existing gap between PBE reaction kinetics and Navier-Stokes equation for the enzymatic hydrolysis of cellulose had been successfully bridged.