Numerical Studies on Cellulose Hydrolysis in Organic–Liquid–Solid Phase Systems with a Liquid Membrane Catalysis Model

The catalytic hydrolysis of cellulose to produce 5-hydroxymethylfurfural (HMF) is a powerful means of biomass resources. The current efficient hydrolysis of cellulose to obtain HMF is dominated by multiphase reaction systems. However, there is still a lack of studies on the synergistic mechanisms and component transport between the various processes of cellulose hydrolysis in a complex multiphase system. In this paper, a liquid membrane catalytic model was developed to simulate the hydrolysis of cellulose and its further reactions, including the adsorption of the liquid membrane on cellulose particles, the consumption of cellulose solid particles, the complex chemical reactions in the liquid membrane, and the transfer of HMF at the phase interface. The simulations show the synergistic effect between cellulose hydrolysis and multiphase mass transfer. We defined an indicator () to characterize the sensitivity of HMF yield to the initial liquid membrane thickness at different reaction stages. decreased gradually when the glucose conversion increased from 0 to 80%, and increased with the thickening of the initial liquid membrane thickness. It was shown that the thickening of the initial liquid membrane thickness promoted the HMF yield under the same glucose conversion. In summary, our results reveal the mechanism of the interaction between multiple physicochemical processes of the cellulose liquid membrane reaction system.


INTRODUCTION
5-Hydroxymethylfurfural (HMF) is a significant chemical raw material that can be obtained from cellulose. 1 However, the physicochemical properties of cellulose are very stable, and there are multitudes of hydrogen bonds within the single chains of the cellulose molecule and between the different chains, leading to cellulose being more difficult to be degraded into monosaccharides. 2 In addition, the preparation of HMF from cellulose is complicated by multitudes of side reactions, and the conversion rate of the raw material is generally low. 3 Therefore, an optimized reaction system is essential for the efficient conversion of cellulose to HMF. 4−6 From the perspective of green chemistry, water is an ideal solvent because it is cheap, non-toxic, and non-flammable. Although cellulose hydrolysis reactions can be effectively carried out using water as the reaction medium, the yield of HMF is generally unsatisfactory due to the instability of HMF in aqueous solution in the presence of an acidic catalyst, which is prone to side reactions. 7 In recent years, researchers have started to use organic solvents as solvents instead of or partially instead of ionic liquids for the preparation of HMF. Polar nonprotonic organic solvents which do not contain hydrogen ions and are less prone to side reactions in the resulting HMF are generally used. In the preparation of 5-HMF from fructose, Benoit et al. 8 partially replaced [BMIM]CL with glycerol as the reaction solvent and obtained a yield of 72% of 5-HMF. However, organic solvents also have some disadvantages, such as difficulty in volatilization and difficult product separation but high boiling points, which is not a general characteristic of all the organic solvents. It appears that neither organic solvents as partial or total substitutes nor water as solvents can give excellent yields of HMF and extract HMF with high purity. Therefore, a biphasic system was proposed to inhibit the redegradation of HMF by combining the aqueous and organic phases, bringing together the reaction to generate HMF and the initial separation of HMF, thus advancing the conversion reaction and allowing for a higher yield of HMF. The group of Dumesic 9−12 at Wisconsin, Madison, focused on the selective variability of HMF when the extracted phase was C3−C6 alcohols, ketones, and furans, making a remarkable contribution to the development of biphasic systems. The results of the research showed that the extraction of HMF was excellent in an organic solvent with carbon atoms of four, and the selectivity of the prepared HMF was the highest when tetrahydrofuran (THF) was used as an extractant.
A reaction system using ionic liquids or high-boiling point organic solvents as reaction media and chromium-containing catalysts showed excellent conversion of HMF from cellulose. 54% yield of HMF was obtained in a combined DMA−LiCl system using HCl and CrCl 3 as catalysts, reported by Binder and Raines 13 Yu et al. 14 catalyzed the degradation of cellulose with CrCl 3 /LiCl bimetallic chloride salts in ionic liquid [EMIM]Cl to obtain 55% HMF. However, the high price of ionic liquids, the high toxicity of chromium chloride, and the high energy consumption for the separation of HMF from the abovementioned systems make it difficult to achieve industrial applications of these studies. Unlike ionic liquids and highboiling point organic solvents, biphasic systems consisting of water and low-boiling point organic solvents are more industrially viable. In the biphasic system, the hydrolysis of cellulose and the formation of HMF take place in the aqueous phase, while the resulting HMF is rapidly transferred to the organic phase to avoid its subsequent degradation. The organic phases commonly used for extraction in biphasic systems are butanol, THF, methyl THF, methyl isobutyl ketone (MIBK), and so forth. However, the yield of HMF obtained in biphasic systems is still relatively low compared to systems consisting of ionic liquids and high-boiling point organic solvents.
Wang et al. 15 obtained only 21% HMF from the conversion of cellulose in a biphasic system consisting of water and butyl phenol, and Yang et al. 16 obtained only 37% yield using cellulose as a feedstock in a biphasic system consisting of water and THF. Shi 17 presented a heterogeneous system for degradation of cellulose to HMF in a biphasic system of water and organic solvents with high concentrations of sulfate. Cellulose hydrolysis occurs on the solid surface. Then, the resulting glucose enters the liquid membrane phase and continues to react to form fructose and HMF. At the same time, HMF diffuses into the bulk phase to reduce the concentration of HMF in the liquid membrane that can significantly inhibit the subsequent side reaction of HMF. The system covers an extensive reaction network including the cellulose hydrolysis, glucose isomerization, HMF dehydration, and complex side reactions. In addition, the system includes three phases and possesses mass transfer at the phase interface. These studies are of great significance to the reaction kinetics and reactive transport phenomena of cellulose to HMF catalyzed by liquid membranes. However, more efforts should be made to gain insights into the mechanisms that influence the physical and chemical phenomena involved in the catalytic hydrolysis of cellulose in liquid membranes to produce HMF. On one hand, as mentioned above, only a few studies have focused on this issue and transport processes are often neglected when designing catalytic systems and reactors. On the other hand, direct measurement and observation of these processes within catalytic liquid membranes are almost impossible.
The lattice Boltzmann method (LBM), which is a mesoscopic numerical algorithm based on a minimal version of Boltzmann's kinetic equation, has become a reliable and efficient simulation technique and received extensive attention. LBM was widely adopted to investigate a variety of physical and chemical phenomena, such as heat transfer, 18,19 phase transition, 20−22 nanofluid, 23−25 and heterogeneous catalysis. 26−28 For multiphase flows, several multiphase LBM models have been developed, of which the Shan−Chen (SC) pseudopotential model is the most widely used due to its superiority over the original concept and its high computational efficiency. Zhang et al. 29 used the LBM to quantify the ionic diffusivity in unsaturated cementitious materials. The results show that the ionic diffusivity was strongly influenced by the degree of water saturation. The simulated relative ion diffusivity as a function of water saturation agreed well with the experimental data obtained from the literature. Model stability and efficiency can be achieved by incorporating a more realistic equation of state (EOS) into the S−C model. 30 A pore-scale model was developed by Chen et al. 31 based on the LBM for the transport of multiphase reactions with phase transitions and dissolution−precipitation processes. The multiphase reaction transport phenomena between liquid and gas phases and the dissolution−precipitation process of salt in a closed envelope were simulated, and the effect of initial envelope size was investigated. Recently, Wei et al. 32 developed a pore-scale multiphase transport model based on the LBM to simulate La to GVL conversion over Ru/C catalysts. Chen et al. 33 developed an SC LB-based pore-scale model by reconstructing the high-resolution porous structure of the cathode catalyst layer to study the reactive transport processes within the catalyst layer nanostructure. These studies demonstrate the power of the LBM in modeling physicochemical phenomena; however, the application of the LBM in heterogeneous catalysis of lignocellulosic biomass has rarely been reported, particularly for complex multicomponent, multiphase catalytic reaction systems such as the hydrolysis of cellulose to produce HMF.
In this context, we aim to develop a numerical model that takes into account multicomponent mass transfer, multiphase flow, phase transitions, and heterogeneous and homogeneous catalytic reactions to enhance the current understanding of the reactive transport phenomena in the production of HMF from cellulose. In this paper, we show a framework for a multiphase reaction transport model based on the SC LB approach, which is used to model the non-homogeneous catalysis of the conversion of cellulose particles to HMF through liquid membrane catalysis in an organic−liquid−solid three-phase system. The multiphase flow is simulated with the SC model with the C−S EOS, and the reactive transport process is treated with the multicomponent mass transport LBM model. A moving boundary condition was considered for the solid− liquid membrane phase interface due to cellulose consumption by the reaction. The effects of initial liquid membrane thickness, reaction temperature, and cellulose particle size were analyzed. We also compare the predictions of the model with experimental results from the literature. 17 Ultimately, effective strategies were proposed to improve the reaction performance and the utilization of the catalytic system based on our conclusion.

MODEL FRAMEWORK
2.1. Multiphase LB Model. In this paper, the Shan−Chen multiphase model (SC model) is used to simulate the liquid membrane reaction and mass transfer process of cellulose particles. The standard LB equation can be expressed as follows where ε f i, is the particle distribution function of component ε along the direction of microscopic velocity e i at site x and time t. τ ε v, is the dimensionless relaxation time related to the kinematic viscosity of the fluid flow. Δt and f i eq are the lattice time step and the equilibrium density distribution function, respectively. The left hand side of eq 1 represents the streaming step, and the right hand side illustrates the collisional relaxation of BGK with the local equilibrium distribution function.
The D m Q n model proposed by Qian is the basic model of the LBM, where m represents the spatial dimension and n represents the number of directions in which the velocity is discrete. In this paper, the two-dimensional D2Q9 model is used, ω i in eq 2 is the weight coefficient, which takes different values in different discrete directions, ω 0 = 4/9, ω − 1 4 = 1/9, ω − 5 8 = 1/36. e i is the lattice discrete velocity, and τ is the single relaxation time and is related to the dynamic viscosity as follows The lattice sound velocity = c c/ 3 s and the lattice velocity = Δ Δ c x t / , where Δt and Δx represent the time step and the space step, respectively.
The fluid density and velocity are obtained from which the zero-and first-order moments of the particle distribution function are calculated.
In the original SC model, Shan and Chen proposed an exponential effective density function, as shown below In order to achieve coexistence between different phase states and adsorption of liquids on solids, two forces are introduced into the LB equation based on the SC model: one is the cohesion between liquid particles on adjacent lattices and the other is the adsorption of a solid surface on a liquid.
The inclusion of a cohesive force that causes phase separation between particles of a particular "species" of fluid, which we call force cohesion, and it is calculated as follows In the equation, F cohe is the interaction force between neighboring fluids, called cohesion and G is the interaction strength constant, whose positive or negative value determines whether the fluid particles attract or repel each other. ψ + x e ( ) i is called the pseudopotential function or effective density, and its value depends on the actual density of the fluid. ω | | e ( ) i 2 is the weight coefficient. For four adjacent grid points , the weight coefficient is 1/3; for four adjacent grid , the weight coefficient is 1/12. When calculating cohesion, we only consider the interaction between adjacent particles. In the D2Q9 model, only the forces of the central particle and neighboring particles in eight directions are calculated.
Based on the results of Yuan and Schaefer, 34 the C−S equation of state is used in this paper, and the Carnahan− Starling equation of state performs better in terms of the twophase density ratio and interface spurious velocity.
The equation of state is used to describe the relationship between fluid pressure, temperature, and density, where = a R T P 0.4963 / 2 c 2 c and = b R T P 0.1873 / c c . In the LB approach, for an ideal gas, since the molecules have no interaction forces, = ρ p 3 . For the SC model, when cohesion within the fluid is considered, the physical quantities are related as follows.
Thus, when we combine eqs 9 with 8 for the C−S state, we can obtain an effective density expression for the variable p and the C−S equation of state.
The adsorption force is calculated as follows In this equation, F adso is the adsorption force of the solid surface to the fluid. 35 G w is the interaction strength constant, which controls the strength of the fluid−solid interaction. + s x e ( ) i is an identification function. It is equal to 1 if the position is a solid phase; otherwise, it is equal to 0, and ρ w is the virtual density of the solid.
By fixing the value of G w and varying ρ w , we can combine the equations for cohesion and adsorption. At the same time, we can treat the solid phase as a liquid phase and set a virtual ACS Omega http://pubs.acs.org/journal/acsodf Article density (which is only used to calculate the adsorption force by interacting with the liquid density).The distinction between the two forces does not need to be taken into account when programming by this method. For the two interacting forces mentioned above, we use the equilibrium velocity correction method proposed by Shan and Chen 36 to incorporate them into the LB equation. In the equilibrium velocity correction method, we consider the influence of force by changing the macroscopic velocity in the equilibrium distribution function. Then, the equilibrium velocity ε u x t ( , eq ) in eq 2 is modified to include the force effect on the ε component as follows In the formula, ε u x t ( , ) eq and * * u x t ( , ) are the equilibrium velocity and common velocity, respectively. The * * u x t ( , ) can be calculated as Substituting eqs 4, 5, and 13 into eq 12, the actual physical velocity is the average velocity before and after the collision and is determined by 2.2. LB Equations for Multicomponent Reactive Transport. To simulate the coupling of passive solute reactive migration with multiphase fluid flow in the LB model, the concentration is modeled using the D2Q9 discrete velocity model, with the evolution equation In the formula, S i is the source term of the chemical reaction, which is related to the chemical reaction rate, τ i,g is the relaxation time related to the diffusion coefficient, and ε g i, eq is the equilibrium function of the concentration distribution and can be calculated by the following formula where ε C is the concentration of εth component. ε J i, and ε i are constants, ε i = 1/2 for 2D simulation, and ε J i, is given by where J i,0 can be selected from 0 to 1 for different diffusion coefficients. The diffusivity and concentration can be obtained by where ε D is the diffusivity of the εth component and c D is a lattice-dependent coefficient, which is equal to 1/2 for the 2D simulation model.

Non-ideal Solute Component
Modeling. Nonideal solutes, such as the HMF studied in this paper, have different solubilities in different solvents. When they are dissolved in multiphase systems, concentration discontinuities occur at the phase interface. Therefore, we need to use methods that enable the solute concentration to pass smoothly through the phase interface.
As shown in Figure 1, the two solvent components σ a and σ b are mixed to form a phase interface, on either side of which are two different phase states, such as the aqueous phase and the organic phase. The whole system contains a non-ideal solute, which is dissolved in both phases. In our simulations, the concentration of the solute transitions smoothly at the phase ACS Omega http://pubs.acs.org/journal/acsodf Article interface, which facilitates the stability of the numerical simulations. In order to achieve a smooth transition of the solute concentration at the phase interface, we introduce a solute−solvent interaction. This interaction affects the distribution of the solute but not that of the solvent. Antoine Riaud 37 uses a multiphase color model to study the diffusion and reaction of dilute solutes in multiphase systems. With reference to the recoloring process of the multiphase color model, we add an additional collision operator to the mass transfer Boltzmann equation to reflect this interaction. An arbitrary function ε δ W X ( ) is chosen to make the solute sensitive to the solvent distribution. The concentration LB equation can be rewritten as In eq 20, the ′ can be written as Combining eqs 20, 22, and, 23, we obtain the following relationship We performed a perturbation analysis on this microdynamic equation and deduced the macroscopic species flux μ is a geometric constant close to 0.150 for D2Q9, and the equation of state of the solute can be derived In eq 24, the vector n is the normal vector δ X of the solvent component distribution. In the pseudopotential model, this can be calculated by the following equation x is a Dirac function; the integral of the force on the diffusion interface leads to the interfacial tension of the sharp interface.
In the collision operator, the value of ε δ W X ( ) is overwhelmingly significant for the solubility of solutes, which we will introduce in detail in specific cases later. γ ε is a physical quantity which determines the degree of solute dissolution in different phases.
2.4. Update of Cellulose Particles. To accurately model the catalytic reactions on the solid surface of cellulose particles, it is a prerequisite to consider the time evolution of the solid phase, especially when the dissolution of cellulose particles involves a mass transfer between the solid and liquid. The amount of dissolved solids is equal to the amount of injected species divided by the difference between the solid concentration and the actual concentration in the cells. In this study, the motion of the solid during fluid flow is not considered, so the volume of the fixed solid is satisfied where ε V , ε V , and ε a are the dimensionless volume, molar volume, and specific surface area of the εth cellulose solids, respectively. In this study, the diffusion of the solute in the solid phase is ignored, and the reaction of cellulose to glucose is considered to occur only at the fluid−solid interface. Each interface node represents a control volume of × 1 1 (lattice unit) (control area in 2D), located in the center of the volume. The initial control volume is the dimensionless volume ε V 0 . According to the equation, the volume is explicitly updated at each time step.
where Δt is the time step and we use the VOP method to track and update the liquid−solid interface. 2.5. Handling Information on Phase-Change Nodes. Figure 2 illustrates the time evolution of the solid−liquid− organic-phase interface. During the interface evolution, the phase of the computational node can be between the liquid and solid phases (due to reactive dissolution of cellulose particles). How to handle the fluid flow and mass transport information associated with the phase transition occurring at the computational node is important to ensure the conservation of mass and momentum in a closed system. For example, the initialization information has various ways to be at the new fluid nodes during the dissolution process and distributed during the multiphase reactive transport, and the information is stored in the newly added entity nodes to the interface. There are three types of nodes in the domain, namely solid nodes, liquid nodes, and organic-phase nodes, and the nodes have dissolution-induced solid to liquid changes. Note that there is no exchange between organic-phase nodes and solid nodes because reactive dissolution occurs only at the liquid−solid interface. When a solid node becomes a fluid node due to dissolution, the density and velocity of that node must be initialized to ensure mass conservation in the closed system and convergence of the simulation. The handling of the flow and concentration information of these nodes is challenging. When the reaction consumes enough cellulose particles to change a solid node into a fluid node, mass conservation, momentum conservation, and species conservation in the system must be guaranteed. Chen et al. 31 propose a phase transition boundary treatment strategy that can guarantee mass conservation and simulation convergence of the system. Specifically, when a solid node becomes a liquid node at time t, the density of this new liquid node (ρ newliq ) is defined as the average density of the nearest neighboring liquid node, and the speed is set as The subscripts "newliq" and "oldsol" denote the new liquidphase node and the old solid-phase node, respectively. As the reaction consumes cellulose to change n nodes in the solid phase to the liquid phase, the liquid-phase density is where σ is a random perturbation to achieve phase transition.
) is the total density of the system at time are new liquid node and old solid node density, respectively. Similar treatment can be adopted for the concentration field.

RESULTS AND DISCUSSION
3.1. Validation of the Model. 3.1.1. Validation of the Multicomponent Multiphase Pseudopotential Model. We use several two-phase flow problems to verify the correctness of the multicomponent multiphase pseudopotential model.
The first problem is a gravity-free suspended droplet in a vapor field, as shown in Figure 3. The number of meshes in the calculation zone is 200 l.u * 200 l.u. All four boundaries are periodic. At the initial moment, a droplet of the radius r 0 is placed in the center of the calculation area and the density of the calculation area is initialized using the following equation  In eq 36: (φ cen = 100, δ cen = 100) is the center of the calculation area. "Tanh" is the hyperbolic tangent function. ρ liq and ρ vap are the theoretical density values at the current temperature (using Maxwell and other density reconstruction methods). In order to verify Laplace's law, we simulate droplets with different radii at temperature = T T 0.72 c . The gas−liquid pressure value is taken from the horizontal middle line AB in (a). Figure 3b shows the pressure distribution on AB when the temperature is = T T 0.72 c and the radius is = r 30 l. u 0 . As shown in (b), the droplet pressure takes the stable pressure value inside the droplet, the vapor pressure takes the stable pressure value outside the droplet, and Δp is the gas−liquid pressure difference. At the gas−liquid interface, the pressure value changes drastically, which is caused by the drastic change in density at the gas−liquid interface. Figure 3c shows Δp at different droplet radii. It can be found that Δp and r 1/ have a linear relationship, which complies with Laplace's law. The slope of the linear fit is 0.0133, which is in good agreement with the theoretical solution of 0.009386. The difference is due to the thermodynamic inconsistency of the LB single-component multiphase model itself.
Another problem in validating the pseudopotential multiphase model is the contact angle. We have created a 200 l.u * 200 l.u grid with a two-dimensional calculation area. The left and right boundaries were periodic boundary conditions, and the top and bottom boundaries were set as solid no-slip boundary conditions. The droplet was placed on the solid wall at the bottom of the computational domain, and the density of the solid was set to the virtual density ρ φ (the virtual density was only used to interact with the fluid density to calculate the adsorption force). Finally, we adjust the virtual density of the solid wall to vary the wettability of the fluid on the solid surface. In the simulation, the virtual density of the solid is defined as follows The coefficient β in the formula determines the wettability of water on a solid surface and takes values from 0 to 1. When β is equal to 0, water is not wettable on a solid surface and when β is 1, water is completely wettable on a solid surface. ρ l and ρ g are the liquid-phase density and the gas-phase density, respectively. When varying the virtual density, the simulation results are shown in Figure 4, where the contact angle gradually decreases from 180 to 0°as the value of β increases from 0 to 1.

Validation of the Mass Transport
Model. In the model, we propose improvements to the mass transfer model for a non-ideal solute concentration distribution in the phase interface region. Therefore, we need to describe and validate the improvements further. In the organic-phase reaction system for the hydrolysis of cellulose to produce HMF, the whole reaction system was divided into three parts, including the solid phase, the liquid membrane phase, and the organic phase. The solutes in the reaction, such as glucose, fructose, and other byproducts, are only dissolved in the liquid membrane. However, there are also substances that dissolve in both the liquid membrane phase and the organic phase, such as HMF, which can dissolve in both phases with different solubilities.
To verify whether the additional collision operator can resolve the interfacial concentration discontinuity, a set of

ACS Omega
http://pubs.acs.org/journal/acsodf Article simulations was carried out based on the SC model, with the center of the simulation domain being the liquid membrane phase (the region where the y-axis is greater than 40 l.u and less than 80 l.u) and the remainder being the bulk phase (the region where the y-axis is less than 40 l.u and greater than 80 l.u). We adjusted the value of ε δ W X ( ) to observe the interfacial concentration distribution of the two solutes. Initially, we set the concentration of glucose in the liquid membrane phase to 1 and the concentration of glucose in the bulk phase to 0. When diffusion reaches equilibrium, glucose is only dissolved in the liquid membrane phase as shown in Figure 5A. In another set of simulations, at the beginning, the concentration of HMF in the liquid membrane was set to 0 and the concentration of HMF in the organic phase was set to 1. When HMF is dissolved in both phases and at diffusion equilibrium, the HMF in the organic phase should diffuse into the liquid membrane     due to the concentration gradient. As shown in Figure 5B, HMF is soluble in both phases at diffusive equilibrium, and the concentration in the organic phase is about five times higher than that in the liquid membrane phase. The solute concentrations in both sets of simulations can be smoothly transitioned in the phase interface region (at y equal to 40 or 60 l.u), where the difference in solubility of the non-ideal solute HMF in the two phases depends on the value of the parameter γ ε in eq 20.

Organic−Liquid−Solid Three-Phase Single-Particle Liquid Membrane Catalytic Model. 3.2.1. Single-Particle Liquid Membrane Catalytic Model.
In the organic− liquid−solid three-phase system, the volume of water added in a general experimental investigation is smaller than the volume of THF. As shown in Figure 6a, the cellulose particles and the aqueous phase containing the catalyst are encapsulated in the organic phase of THF. Khazraji et al. 38 suggest that due to a large number of hydroxyl groups in the cellulose molecule, it has strong hydrophilic properties, which means that it is also highly hygroscopic in multiphase solvent systems. At room temperature, the surface of cellulose is also covered with a layer of water molecules. The polarity of the water molecules is    The cellulose particles will inevitably preferentially adsorb the water phase and form an acidic liquid membrane on the surface of the cellulose particles. We can assume that the surface of the cellulose particles is completely covered by a water membrane; this liquid membrane containing a high concentration of the catalyst can effectively catalyze the cellulose hydrolysis reaction to prepare HMF, while the HMF molecules produced by the reaction are rapidly transferred to the organic phase to avoid subsequent side reactions. In this three-phase system, the organic phase has three functions. The first is to act as a solvent for extracting HMF and transferring it out of the reaction phase; the second is to act as a dispersant for the cellulose particles and the aqueous phase; the third function is a storage area for HMF. In addition, water also has three functions: one is to act as a green solvent, dissolving the catalyst and providing the reaction environment required for the reaction; another is to act as a reactant, hydrolyzing the cellulose to produce glucose, and the last is to change the volume ratio of the aqueous phase to the organic phase.
The hydrolysis of cellulose to HMF in the organic−liquid− solid system is a complex reaction. First, the hydrolysis of cellulose occurs on the surface of the particles, producing glucose, and then, the resulting glucose will diffuse into the liquid membrane for a subsequent series of reactions. The conversion of glucose to HMF can be divided into two steps. In the first step, glucose is converted to the more chemically active fructose by isomerization. Although there is still some debate about the reaction mechanism of the isomerization process, it is widely accepted that fructose is used as an intermediate product in the reaction process. 39−46 The second step is the removal of three molecules of water from the fructose to produce HMF, which is chemically active in acidic aqueous solutions and is prone to a variety of side reactions. Therefore, timely transfer of HMF and inhibition of HMF side reactions in the liquid membrane are key aspects to improve

ACS Omega
http://pubs.acs.org/journal/acsodf Article HMF yield and achieve efficient preparation of HMF from cellulose. In liquid membranes, not only HMF undergoes side reactions but also the intermediate products glucose and fructose are accompanied by a correspondingly large number of side reactions, producing humins, levulinic acid, and formic acid. 39−46 In this reaction, the catalyst comes into contact with the cellulose very easily, so the reaction system combines the advantages of both homogeneous and heterogeneous catalysts. Cellulose in liquid membrane catalysis itself acts as a catalyst carrier (Table 1).
Referring to the reaction network in Figure 6b   where k a and − k a are the rate constants for the isomerization of glucose to fructose and fructose to glucose, respectively; * k a is the net rate constant for isomerization ( * = − − k k k a a a ); k d is the rate constant for the dehydration of fructose to HMF; k g is the rate constant for the rehydration of HMF to FA/LA; k b , k e , and k h are the rate constants for the degradation of glucose, fructose, and HMF to FA and other products (oligomers, humic substances, etc.), respectively; and k c , k f , and k i are the rate constants for the polymerization of glucose, fructose, and HMF to form humin substances, respectively. The values of the reaction rate constants are quoted from Tang 47 and Yan. 48 We can obtain the kinetic constants for the reactions at different temperatures through the Arrhenius equation  In the equation, R is the molar gas constant, T is the reaction temperature, E a is the apparent activation energy, and A is the pre-finger factor. The values of the reaction rate constants at different temperatures are shown in Table 2. The following equations are used to connect the actual physical parameters to the dimensionless physical parameters, the particle diameter is 4.2 × 10 −4 m, and the kinematic viscosity ϑ is 1.006 × 10 −6 m 2 s −1 .
We investigated the organic−liquid−solid three-phase system for the hydrolysis of cellulose to produce HMF using a previously validated numerical model. To simplify the numerical simulation work, a liquid membrane catalytic model was developed under the assumption that the cellulose particles are rigid spheres and that the cellulose particles will continue to dissolve during the reaction, as shown in Figure 7. The simulation area is 250 l.u*250 l.u, which is a three-phase system. Initially, we placed spherical cellulose particles with a particle size of 50 l.u (orange area) in the center of the calculation region. Due to the hydrophilic nature of the cellulose particles, a liquid membrane (pale-yellow area) with a thickness of 10 l.u is adsorbed on the surface of the cellulose particles. The outermost layer is the organic phase (light-blue area). In the middle of the pale-yellow and light-blue areas is the phase interface between the liquid and organic phases. Periodic boundary conditions were used for the flow and concentration fields at all four boundaries of the calculation zone, and half-way bounce-back boundary conditions were used for the flow and concentration fields at the solid surface. The simulation begins with the hydrolysis of cellulose at the particle surface to produce glucose, and the cellulose particles continue to dissolve and shrink. The glucose produced then diffuses into the liquid membrane and participates in subsequent reactions to produce fructose, HMF, and other byproducts. The HMF produced during this process will diffuse into the organic phase, which reduces the concentration of HMF in the liquid membrane and avoids its subsequent side reactions.
3.2.2. Results of the Model. In the reaction of HMF formation from the hydrolysis of cellulose particles in the organic−liquid−solid-phase system, we focused on the variation of the three products, glucose, fructose, and HMF. Figure 8 shows a cloud chart of the distribution of the concentrations of HMF, glucose, and fructose throughout the reaction system with reaction time. At t = 10,000 t.s, only the HMF concentration near the organic−liquid membrane phase interface is higher in the entire organic phase; however, the HMF concentration is lower at the boundary of the calculated region. In addition, the HMF concentration in the liquid membrane phase is lower than the HMF concentration near the phase interface in the organic phase. The concentration of HMF in the organic phase spreads outward with the phase interface and shows a certain concentration gradient, with a relatively high concentration of HMF nearer to the phase interface. Glucose and fructose are distributed in the liquid membrane phase, where the concentration is higher, and the distribution is uniform and no concentration gradient occurs. The concentration of fructose and glucose decreases as the reaction proceeds, with fructose being more uniformly distributed in the liquid membrane phase and glucose being more highly distributed near the surface of the cellulose particles, where a concentration gradient has occurred. The concentration of glucose is lower near the interface between the organic and liquid membrane phases. As the reaction proceeded to t = 30,000 t.s, the concentration of HMF in the organic phase was higher near the phase interface and also increased at the edge of the organic phase. When the reaction proceeds to t = 40,000 t.s, the concentration of HMF in the organic phase is higher and the concentration at the boundary of the calculated region is already higher than that in the liquid membrane. When the reaction proceeds to t = 50,000 t.s, the concentration of HMF in the whole organic phase is already high and the low concentration of HMF in the liquid membrane can inhibit the occurrence of HMF side reactions, which is conducive to the increase in HMF yield.
From the cloud chart in Figure 8, we can find the product concentration distribution characteristics of glucose, fructose, and HMF when they proceed through the entire reaction system. Glucose and fructose are distributed in the liquid membrane phase, where the glucose concentration is relatively low, further away from the solid−liquid phase interface. The glucose concentration distribution shows a linear relationship with the distance from the solid−liquid phase interface. The fructose concentration basically maintained the uniform distribution pattern in the liquid membrane phase. We investigate their concentration distribution comprehensively. In particular, glucose and fructose are only distributed in the liquid membrane phase, whereas HMF, a non-ideal solute, is distributed in both the liquid membrane phase and the organic phase. Figure 9 illustrates the normalized concentration distributions of glucose, fructose, and HMF during the reaction. Since the main object is the concentration distribution within the liquid membrane, the normalized concentration is defined as = C C C / N glucosemax . In this equation, C glucosemax represents the maximum value of the glucose concentration throughout the reaction. We intercepted the reaction concentration curve on a cross-section of y = 125 l.u to investigate the distribution of the concentrations of the three products in the x-direction. The chart is divided into three parts, namely the solid phase, the liquid membrane phase, and the organic phase. Since we assume that the cellulose particles are rigid spheres, the concentration of HMF in the solid phase is zero and does not need to be discussed.
From Figure 9, we can investigate the distribution characteristics of glucose, fructose, and HMF in the whole reaction system. Glucose and fructose were only distributed in the liquid membrane phase, and the concentration of glucose was lower further away from the cellulose particles and showed a certain linear relationship with the distance from the phase interface. The concentration of fructose remained essentially uniform in the liquid membrane phase and decreased as the reaction time progressed. When the reaction proceeded up to t = 10,000 t.s, the concentration values of HMF in the liquid membrane phase or organic phase were low, and at t = 30,000 t.s when the reaction proceeded, the concentration of HMF started to increase in the liquid membrane phase, and the concentration of HMF was higher in the region close to the cellulose particles. At this point, because of the transport of HMF to the phase interface, concentration diffusion causes a lower concentration at the phase interface, while at the organic−liquid membrane phase interface, the concentration of HMF shifts dramatically because of the interfacial transport mechanism of the non-ideal solute. We observed that at the boundary of the calculated zone, the concentration values of HMF increased with reaction time and the high concentration of HMF in the organic phase favored the increase in HMF yield. Crucially, the low concentration of HMF in the liquid membrane phase significantly inhibited the subsequent side reactions of HMF and increased the yield of HMF significantly. At t = 50,000 t.s, the cellulose particles gradually decreased in size as the reaction proceeded further, and the peak HMF concentration in the liquid membrane phase of the reaction system did not increase. As the cellulose particles shrank further, the liquid membrane phase of the overall reaction system continued to increase. In the liquid membrane phase, the concentration of HMF is relatively low further away from the cellulose particles.
3.3. Mechanism of the Initial Liquid Membrane Thickness in the Reaction System. The initial liquid membrane thickness ξ plays an important role in the reaction system as a measure of the liquid membrane phase. By analyzing the simulation data, we obtained relevant conclusions about the role of the initial liquid membrane thickness in influencing the mechanism of the reaction system.
3.3.1. Conversion of Glucose. With different initial liquid membrane thicknesses, there is some variation in the conversion of glucose throughout the reaction system. As can be seen from Figure 10, when the reaction proceeds, the shortest reaction time is required for the conversion of glucose to reach 100% when the initial liquid membrane thickness is ξ = 3 l.u. When the initial liquid membrane thickness ξ continues to increase, the time required for the conversion of glucose to reach 100% also continues to increase. Because when the initial liquid membrane thickness is thinner, the concentration of catalyst H + in the liquid membrane is higher. Therefore, the reaction rate of glucose is relatively fast, leading to a more rapid increase in the conversion of glucose.
3.3.2. Sensitivity of HMF Yield to Initial Liquid Membrane Thickness. After analyzing the mechanism of the effect of different initial liquid membrane thicknesses on the conversion of glucose, we proceeded to investigate the variation curves of the yield of HMF with the glucose conversion proceeding at different initial liquid membrane thicknesses for a cellulose particle size of 50 l.u. Here, we propose Ω a to represent the indicator for the sensitivity analysis of HMF yield to initial liquid membrane thickness. We take the HMF yield at each glucose conversion with an initial liquid membrane thickness of ξ = 3 l.u as a benchmark and establish the relative increment of HMF yield for other initial liquid membrane thicknesses relative to 3 l.u. The equation for Ω a is as follows Based on the data related to Ω a , we made the chart in (a) of Figure 11. According to the trend of the chart (a), we divide it into three stages, 1, 2, and 3. Figure 11b shows the trend of the Ω a change in stage 1. In the first stage, the value of Ω a gradually decreased as the conversion of glucose increased from 0 to 80%. Furthermore, the value of Ω a increases as the initial liquid membrane thickness ξ continues to increase. This indicates that in stage 1, the larger the initial liquid membrane thickness is, the higher the value of Ω a is relatively. The thickening of initial liquid membrane thickness plays a facilitating role in the growth of HFM yield. It shows that the thicker the initial liquid membrane thickness is at this stage, the more sensitive the effect on the HFM yield is because at the beginning of the reaction, larger initial liquid membrane thickness in the reaction leads to a relatively larger volume of the liquid membrane phase, which also means that the reaction area for the generation of HMF from glucose in the liquid membrane phase is relatively larger. The growth rate of the yield of HMF also increases with the increase in the initial liquid membrane thickness. In stage 2, when the conversion of glucose increased from 80 to 96%, we could obtain that Ω a was below 0% at this time, indicating that the highest yield was achieved at the initial liquid membrane thickness of 3 l.u in the second phase and Ω a decreased more rapidly with the increase in the initial liquid membrane thickness. It shows that the initial liquid membrane thickness played a suppressive role in the yield of HMF at this stage as the initial liquid membrane thickness increased. This is because when the reaction proceeds for a period of time, the catalyst H + concentration in the liquid membrane phase decreases with the increase in the initial liquid membrane thickness and the reaction rate decreases, which leads to a decrease in the yield of HMF. Entering the third stage, the glucose conversion was between 96 and 100% at this time. When the glucose conversion was close to 100%, the value of Ω a showed an increasing trend at this time, which was due to the fact that the side reactions occurred more significantly after the complete ACS Omega http://pubs.acs.org/journal/acsodf Article reaction of glucose when the initial liquid membrane thickness was thinner, leading to a decrease in HMF yield. As a result, we derived a mechanism for the regulation of HMF yield and the reaction system at different reaction stages by different initial liquid membrane thicknesses. We can obtain relatively high HMF yields at different reaction stages by varying different initial liquid membrane thicknesses. 3.3.3. Selectivity of HMF. After analyzing the mechanism of the initial membrane thickness on the HMF yield at different reaction stages, we then analyzed the effect on the HMF selectivity. Figure 12 shows the curve fitted to the HMF selectivity as a function of glucose conversion from 70 to 95%. It can be found that the HMF selectivity tends to increase and then decrease as the conversion rate of glucose increases. It is found that the HMF selectivity increases as the initial membrane thickness increases from 3 to 7 l.u, and then, the HMF selectivity reaches a maximum at an initial membrane thickness of 7 l.u.Then, the selectivity of HMF gradually decreases as the initial membrane thickness increases from 7 to 15 l.u. At this point, we can say that there exists an optimal liquid membrane thickness of 7 l.u that keeps the selectivity of HMF at its maximum.
3.4. Role of Cellulose Particle Size in the Reaction System. After we investigated the mechanism of the influence of different initial thicknesses in the overall reaction system in the previous section, we proceeded to investigate the role of different cellulose particle sizes r in the reaction. As shown in Figure 13, the variation curves of glucose conversion, HMF yield, and fructose selectivity at different cellulose particle sizes are shown. In Figure 13a, the reaction time required to reach 100% glucose conversion will be prolonged when the cellulose particle size keeps increasing because the increase in cellulose particle size corresponds to the decrease in H + particle concentration, leading to reaction rate decreases, and the glucose conversion rate decreases subsequently. The maximum value of HMF yield decreases continuously as the cellulose particle size increases, and the yield of HMF at each moment decreases with the increase in cellulose particle size in Figure  13b. The selectivity of fructose was also correlated with the cellulose particle size in Figure 13c, and the selectivity of fructose was close to 0% when the reaction proceeded to t = 50 t.s. It was found that the selectivity of fructose decreased faster with the increase in cellulose particle size.
3.5. Effect of Reaction Temperature in the Reaction System. After investigating the roles played by the initial liquid membrane thickness and cellulose particle size in the reaction system, we proceeded to investigate the effects of different reaction temperatures on the reaction system. Figure  14 shows the cloud chart of the distribution of the three main product concentrations at different reaction temperatures at the reaction time of t = 30,000 t.s. When the temperature was increased from 413 to 433 K, the concentration of HMF increased significantly and was concentrated at the organic− liquid phase interface, while the concentrations of fructose and glucose decreased with the increase in temperature, and the distribution was more uniform at this time. When the temperature was increased from 433 to 453 K, the distribution of HMF in the organic phase was more concentrated and the concentration continued to increase. The concentrations of fructose and glucose continued to decrease and a concentration gradient appeared, with a relatively high concentration near the solid−liquid phase interface.
We observed changes in the concentration distribution of the three main products by varying different reaction temperatures, which similarly affect the conversion of glucose and the yield of HMF. As shown in Figure 15, Figure 15a shows the curves of glucose conversion with time for different reaction temperatures. When the reaction temperature is higher, the glucose conversion grows faster and the time required to reach the maximum conversion is shorter, while at lower reaction temperatures, the glucose conversion grows slower and the time required to reach the maximum conversion is substantially longer, and the conversion of glucose cannot reach 100% at lower temperatures. Figure 15b shows the curves of HMF yield with reaction time at different reaction temperatures. An optimal reaction temperature of 453 K was exerted to get the highest yield of HMF in Figure 15b. After the reaction time exceeded t = 20 t.s, we can keep the HMF yield relatively high by regulating the reaction temperature.
We compare the simulation results with the experimental results, 17 and basically, the trend of the change is in agreement, as shown in Figure 16. In the experimental results, only 35.8% of HMF was generated at the lower temperature of 413 K, and when the temperature was increased, the yield of HMF increased significantly. When the temperature was increased from 413 to 433 K, the yield of HFM also increased to 53.2%. Then, when the temperature was increased to 453 K, the HMF yield decreased to 45%. This is also in line with the trend of our simulated results.

CONCLUSIONS
An LBM-based organic−liquid−solid multiphase mixing model was developed to simulate physicochemical problems involving organic−liquid−solid-phase states and to model the liquid membrane catalytic system for the reaction of solid cellulose particles to generate HMF. The multiphase liquid membrane catalytic reaction model can capture well the multicomponent mass transfer, homogeneous and heterogeneous reactions, and reactive dissolution processes of the solid phase. The effects of initial liquid membrane thickness, reaction temperature, and cellulose particle size on the multicomponent multiphase reaction transport were investigated and discussed in detail. The following conclusions were drawn: I. When the reaction temperature was increased from 413 to 453 K, the presence of an optimum reaction temperature of 433 K allowed the yield of HMF to remain at the highest value for a longer reaction time. II. Under the conditions of constant reaction temperature and initial liquid membrane thickness, the yield of HMF gradually increased as the cellulose particle size increased from 40 to 60 l.u; however, the reaction time required to reach the maximum glucose conversion gradually increased with the increase in cellulose particle size. III. At different initial liquid membrane thicknesses, as the reaction proceeded, the shortest reaction time was required to achieve 100% conversion of glucose when the initial liquid membrane thickness was thin. An indicator Ω ( ) a was defined to characterize the sensitivity of HMF yield to the initial liquid membrane thickness at different reaction stages. Ω a decreased gradually when the glucose conversion increased from 0 to 80%, and Ω a increased with the thickening of the initial liquid ACS Omega http://pubs.acs.org/journal/acsodf Article membrane thickness. It was shown that the thickening of the initial liquid membrane thickness promoted the HMF yield under the same glucose conversion. When the glucose conversion rate increased from 80 to 96%, Ω a decreased with the increase in the initial liquid membrane thickness, and the results indicated that the thickening of the initial liquid membrane thickness inhibited the increase in HMF yield at the same glucose conversion. The selectivity of HMF tends to increase with the increase in glucose conversion and then decreases, and there is an optimal membrane thickness of 7 l.u to keep the selectivity of HMF at the maximum.