Article pubs.acs.org/IECR
Modeling of Mass Transfer and Reactions in Anisotropic Biomass Particles with Reduced Computational Load Zheng Liu,*,† Ville Suntio,† Susanna Kuitunen,† Jonas Roininen,‡ and Ville Alopaeus† †
Department of Biotechnology and Chemical Technology, School of Chem. Tech, Aalto University, P.O.B. 16100, FI-00076 Espoo, Finland ‡ Neste Jacobs Oy, P.O. Box 310, FI-06101 Porvoo, Finland ABSTRACT: Accurate description of mass transfer and chemical reactions in anisotropic biomass particles is crucial when formulating truly predictive biorefinery process models. Due to the natural anisotropy of biomass particles, these can be accurately modeled with partial differential equations only by considering the mass transfer in all three dimensions. However, a large amount of memory and CPU time is needed to solve this three dimensional (3D) model. In order to reduce the computational load, this article introduces a method to simplify the 3D model into a 1D model. Furthermore, the traditionally used low order finite difference method is replaced by the high order moment based weighted residual method. Using the simplified 1D modeling approach and the new numerical method, a much lower number of variables is needed than used when solving the 3D model with the finite difference method, meanwhile providing similar average concentration profiles as the 3D model.
1. INTRODUCTION Biomass, representing an abundant carbon-neutral resource, includes wood residues, agricultural wastes, algal, etc. During the first half of 20th century, many industrial materials such as dyes, solvents, and synthetic fibers were already made from biomass such as trees and agricultural crops. After that, many of these biobased chemical products had been displaced by petroleum derivatives.1 Due to the climate change and growing demand on limited petroleum resources, the interest on biomass as the raw material has been restored steadily during this century especially for the production of bioenergy, liquid transportation fuels, and a wide variety of precursor chemicals. The manufacturing concept for converting biomass to various valuable products is referred as biorefinery. For example from wood and other lignocellulosic biomass, liquid fuels such as bioethanol, biodegradable plastics such as polyhydroxyalkanoates, and acetic acid can be produced by biorefinery processes.2 For processes handling, either inorganic catalyst particles or biomass particles, both mass transfer and reactions are critical and relevant phenomena. For design and optimization of novel biorefinery processes utilizing biomass particles, phenomenon based modeling is a very efficient tool. However, due to the anisotropic nature of biomass particles and complex reactions taking place inside the particles, the resulting models usually include a high number of variables consequently leading to long simulation times. The aim of this work is to reduce the computational burden of solving the concentration profiles developing inside the biomass particles. The kraft pulping of wood chips was selected as the study object of this work, to represent various biomass particles and related biorefinery processes. This is because the kinetic models of wood chips are already relatively well studied and wood chips are one of the most common and widely used biomass particles in different fields. Also for the reason that like © 2014 American Chemical Society
the modeling of wood chip pulping, modeling of any other biorefinery processes utilizing biomass particles may encounter similar challenges related to shape, anisotropy, and number of components involved in reactions. In general, kraft pulping of wood chips is to solubilize lignin by breaking up its ether bonds and at the same time minimize the degradation of carbohydrates. During the pulping process, the main reactions between lignin, carbohydrates, and the inorganic reactants take place inside the porous structure of the solid matrix of wood chips. Since it is generally known that chemical reaction rates are faster than diffusion rates in wood chips,3,4 mass transfer of the inorganic reactants plays an important role in the overall rate of impregnation and pulping processes. Various mass transfer oriented models of pulping processes have been proposed and implemented.5−8 Among them, the models taking into account the mass transfer in all three dimensions (3Ds) of rectangular anisotropic wood chips are considered as the most accurate mathematical description of the process. Similar to other biomass particles, due to the anisotropic nature of wood chips and a large number of components and reactions involved, modeling the development of concentration profiles inside the wood chips with the 3D model is highly complex and the computational burden is heavy. In this paper, the terms geometrical shape parameter and particle characteristic dimension, which account for anisotropy and dimensions of wood chips, are first introduced. By these means, the rigorous 3D wood chip model is simplified into a one-dimensional (1D) model. Accuracy of the modeling is maintained while the number of variables to solve the model dramatically decreases. Traditionally low-order finite difference Received: Revised: Accepted: Published: 4096
October 10, 2013 February 6, 2014 February 18, 2014 February 18, 2014 dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research
Article
⎡ d2C dCi d2Ci d2Ci ⎤ i ′ ′ ⎥ + Sr E E = Di⎢Ex′ + + y z dt dy 2 dz 2 ⎦ ⎣ dx 2
method is used for solving these kinds of partial differential equations. A high order schememoment based weighted residual method (moment method)was applied recently to solve dynamic reactor models and chromatographic models.9−11 In the present paper, the implementation and general formulation of the moment method are introduced to solve the 1D model. By this high order method, further decrease in the number of variables and computational burden is expected. The validation of the 1D model, solved either with finite difference or moment method, is done by comparing the average concentration profiles inside the chip to those obtained with the 3D model, which is considered as the reference case in this study.
(1)
Because of symmetry, the concentration gradients vanish at the center of the chip. The flux at the surface of the chip is considered to be continuous. This means that the flux across the film and the surface of the chip are equal to each other. The boundary conditions at the center (eq 1a) and at the surface of the chip (eq 1b) are the following: x = y = z = 0;
2. RESEARCH METHOD In this work, the mass transfer and reactions taking place inside the anisotropic biomass particles are modeled with three approaches: (1) The 3D model considers the anisotropic feature of the rectangular particle and accounts for different mass transfer resistance in all three dimensions. (2) In the homogeneously mixed chip phase model, the concentration profiles inside the particle are ignored, i.e., the phase inside the particle is homogenously mixed and all the mass transfer resistances are lumped into the film on the particle surface. (3) The 1D model is a dimensional reduction of the 3D model. It is solved with either finite difference or moment method to simulate the concentration profiles inside the particle. Compared to other models, the 3D model is the most rigorous approach to describe the mass transfer phenomena in anisotropic particles. Average concentration profiles obtained with the 1D model and the homogeneously mixed chip phase model are validated by comparing them to those obtained with the 3D model. 2.1. 3D Model. In this work, the 3D model is formulated for wood chips due to their abundance as a raw material and importance in variety of applications related to biorefinery processing. However, it should be noted that the present approach is by no means limited to wood as the starting material. A number of fundamental assumptions to construct the model include the following: The wood chips are porous and anisotropic. The reaction occurs only in the wood chip phase. The mass transfer in wood chips happens via diffusion. The model includes the internal pore diffusion resistance and the external film diffusion resistance. The wood chips ordinarily undergo structural changes, e.g., in the delignification process. In this model, the structure of the wood chips is assumed to be unchanged because its structure change process is slow. Consequently the parameters, representing the chip structure, i.e., porosity and tortuosity, remain constant. The detailed derivation and assumptions of the model were described thoroughly in the article of Grénman.8 The distinction between the present model and Grenman’s model is in the definition of boundary conditions. The external film mass transfer resistance coefficient “k” is used in the present model. In Grenman’s model, the surface concentrations of the particle phase are considered to be the same as of the bulk phase; therefore, the external film diffusion resistance was neglected in their study. The model of the wood chip phase is
dCi dCi dCi = = =0 dx dy dz
x = X,
kiτyz ∂Ci = (C bi − Ci , x = X ) ∂x εyzDi
y = Y,
kτ ∂Ci = i xz (C bi − Ci , y = Y ) ∂y εxzDi
z = Z,
kiτxy ∂Ci = (C bi − Ci , z = Z) ∂z εxyDi
where E′ is the effective capillary cross-sectional area as εyz εxy ε Ex′ = , E′y = xz , Ez′ = εpτyz εpτxz εpτxy
(1a)
(1b)
(1c)
Sr is the reaction term. The simple reaction involved in this study is: C + W ⇒ P with the reaction rate law Sr = KrCcCw. C represents the inorganic reactant, e.g. alkali; W represents the lignin existing in the wood chips; P represents the product of the reaction, e.g. the dissolved lignin. This reaction is very simplified representation of the delignification in pulping, because the main focus of this study is to simplify the 3D model and to develop the numerical strategies to solve the model. In the following discussion concerning other modeling approaches, the reaction term Sr remains identical. The mass balance of the liquid bulk phase is d C bi (1 − εb) = [Gyz(Nix) + Gxz(Niy) + Gxy(Niz)] dt εbVp
(2)
where εb is the bulk phase volume proportion in the bulkparticle phase system. VP is the wood chip volume. G stands for integration over the whole surface area, e.g., Gyz(Nix) = ∫ Y0 ∫ Z0 Nix dy dz. Ni are the fluxes in different directions at the surface of the chip: εyz Nix = − ki(C bi − Ci)x = X τyz Niy = − Niz = −
εxz ki(C bi − Ci)y = Y τxz εxy τxy
ki(C bi − Ci)z = Z (2a)
2.2. Homogeneously Mixed Chip Phase Model (HMC Model). Except the assumptions in the 3D model, one important additional assumption here is that the concentration profiles inside the chips are ignored and the chip phase is considered to be homogeneously mixed in this most simplified modeling approach. This method is often used in the modeling of digesters.7 The mass transfer resistance exists only at the 4097
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research
Article
these parameters is introduced in the following section. The liquid bulk phase mass balance equation is
outer surface of the chip. In this case, the previous model eq 1 and 2 are simplified as eq 3 and 4 presented below. The mass balance of the wood chip phase is
d C bi (1 − εb) = (NA i α ), dt εbVα εp where Ni = − ki(C bi − Ci , R = L) τ
dCi 1 = [ki , xyεxy(C bi − Ci)Axy + ki , xzεxz(C bi − Ci)Axz dt εpVp + ki , yzεyz(C bi − Ci)A yz ] + Sr
(3)
2.3.1. Estimation of the 1D Model Parameters. One goal in this study is to test whether the 1D model as eq 6 gives accurate enough description of the system, when the geometrical shape parameter is estimated for a rectangular wood chip simultaneously accounting for anisotropy of the wood. This idea is based on the work of Neretnieks.13 In the present model, it is assumed that E′ remains constant throughout the simulation. Equation 1 (leaving out the reaction term Sr) is divided by Di· Ex′:
The mass balance of the liquid bulk phase is d C bi (1 − εb) =− [ki , xyεxy(C bi − Ci)Axy dt εbVp + ki , xzεxz(C bi − Ci)Axz + ki , yzεyz(C bi − Ci)A yz ] (4)
The parameters kxy, kxz, and kyy represent the lumped mass transfer resistance coefficients in three directions of the wood chip. According to the following empirical correlations,5 the lumped values of these coefficients can be calculated from the diffusion coefficient D and the external film diffusion coefficient k can be calculated as defined in the 3D model. 1 ki , xy 1 ki , xz 1 ki , yz
=
Lτxy 1 + 3Di ki
=
Lτ 1 + xz 3Di ki
=
⎡ d2C E′y d2Ci Ez′ d2Ci ⎤ 1 dCi ⎥ = ⎢ 2i + + Ex′ dy 2 Ex′ dz 2 ⎦ Di ·Ex′ dt ⎣ dx
(5)
D 1 ∂ ⎛ Ω ∂Ci ⎞ ∂Ci ⎜R ⎟ + Sr = i Ω ∂t τ R ∂R ⎝ ∂R ⎠
Y′ =
Z′ =
(10)
R = 0,
Ω=
(11)
Z Ez′ /Ex′
(12)
Aα L−1 Vα
(13)
Aα is the computational surface area of the particle, estimated with eq 14. Aα = 2·(XY ′ + XZ′ + Y ′Z′)
(14)
Vα is the computational volume of the particle, estimated in eq 15.
Vα = XY ′Z′
(15)
The characteristic dimension, L, equals to the smallest dimension among X, Y, and Z divided by two. Tortuosity τ of eq 6 equals to the smallest tortuosity value among τxy, τxz, and τyz. The general methodology presented in this paper can be also extended to irregular wood chips or biomass particles, for which the geometrical shape parameter, Ω, is calculated as in the article of Burghardt and Kubaczka.12 2.3.2. Moment Method to Solve the 1D Model. Before the implementation of the moment method is presented, the model in eq 6 and the initial and boundary conditions in eq 6a are transformed to dimensionless form by introducing the following dimensionless terms:
(6)
∂Ci = 0; ∂R
kτ ∂Ci = i (C bi − Ci , R = L) εpDi ∂R
Y E′y /Ex′
These new dimensions (X, Y′, and Z′) are used for the computation of geometrical shape parameter, Ω, which is defined as eq 13:12
with the initial and boundary conditions
R = L,
z Ez′ /Ex′
The new dimensions of the chip are then:
2.3. 1D Model. Although the 3D model describes the phenomena taking place inside the chips accurately, its computational burden is however exceedingly heavy. In order to lower the memory and CPU demand, a 1D model is developed. The starting point of the 1D model is the generalized PDE for particles as in eq 6, for which mass transfer needs to be considered only in one dimension. The parameter Ω is the geometrical shape parameter, which has been used in catalyst field to generalize the effective factor of any shape of a catalyst pellet.12 The assumption of constant porosity and tortuosity presented in section 2.1 is critical because it is the foundation to estimate the parameter Ω as introduced in section 2.3.1. In this study, we intend to generalize the rectangular shape of the wood chip and other properties related to mass transfer (e.g., the anistropic tortuosity and porosity). With proper choice for Ω, the equation can be used for slab (Ω = 0), cylinder (Ω = 1), and sphere (Ω = 2).
t = 0, C = [CCC PC W ];
(8)
The spatial variables are then changed by doing the following substitutions: y y′ = E′y /Ex′ (9)
z′ =
Lτyz 1 + 3Di ki
(7)
(6a)
Cb is concentration in the liquid bulk phase. The parameter L is the characteristic dimension of the particle. The estimation of 4098
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research ci =
Article
Ci tD Di k Lτ R , θ = ref , r = , ηi = , Bii = i , 2 Cref L Dref τ εPDi L
kr =
In order to investigate the feasibility of the 1D model and the moment method, we compare the average concentration profiles in the particle phase and in the liquid bulk phase from the 1D model and the 3D model. One advantage of the moment method is that the average concentration of the particle phase can be conveniently calculated from the Ωth moment. Next the relationships between the average concentration of the particle phase and the liquid bulk phase with the Ωth moment are presented. The Ωth moment for the time derivative of the average concentration is written as
K rCref L2 Dref
where Cref = 1.0 mol/L; Dref is the diffusion coefficient value of component C. The dimensionless formulation of the model in eq 6 is then ∂ci 1 ∂ ⎛ Ω ∂ci ⎞ ⎟ + k rccc w ⎜r = ηi Ω ∂θ r ∂r ⎝ ∂r ⎠
dm Ω = dθ
(16)
with the dimensionless initial and boundary conditions θ = 0, c = [ccc pc w ]; r = 1,
dθ
=
∫0
1
(16a)
∂c j r dr ∂θ
dθ
=
∫0
η
dθ
=η
∫0
1
⎛ ∂c ⎞ r j −Ω d⎜r Ω ⎟+ ⎝ ∂r ⎠
dθ
= ηBi(cb − cr = 1) − +
∫0
1
r Ω dr =
∂cave 1 (Ω + 1) ∂θ
∫0
1
k rccc wr j dr
∫0
1
k rccc wr j dr
∫0
1
η(j − Ω)r j − 1
k rccc wr j dr
(24)
εp(1 − εb) dcave εp(1 − εb) dc b dm Ω =− =− (Ω + 1) εb εb dθ dθ dθ (25)
Equation 25 is the general expression of the bulk phase for a random Ω value. If the value of Ω is set as 2 and the reaction term in eq 21 is ignored, the bulk phase model becomes identical with the classical chromatographic general rate model14 without the convection term. Other implementation details of the moment method can be found in the literature.9−11 2.4. Quantitative Evaluation of Modeled Profiles. In this work, the concentration profiles simulated with different models and numerical schemes are compared to the 3-D model solved with the finite difference method. Except the visual difference, the average quantitative deviations for all the simulation points of these profiles from the 3D model are also calculated and presented as
(18)
(19)
∂c dr ∂r
1 dev = 2mn (20)
∫0
1
k rccc wr Ω dr
n
m
∑∑ i=1 j=1
(Yi , j − Yi , j ,3D)2 Yi , j ,3D2
(26)
where, 2 represents number of phases, i.e. the particle phase and bulk phase; m is total number of simulation points; n is number of components; Yi,j is the jth simulated average concentration of component i by various models and numerical schemes; Yi,j,3D the jth simulated average concentration of component i by the 3D model The relative errors of individual simulated concentrations by various models and numerical schemes compared to the 3D model are obtained from
From eq 20, only when j = Ω, the second term on the righthand side vanishes. The equation of Ωth moment appears as dm Ω = ηBi(cb − cr = 1) + dθ
(23)
The bulk phase concentration is related to the average concentration of the particle phase and Ωth moment as in eq 25.
Integrating by parts of the mass transfer term and inserting the boundary conditions, the final expression of the model moment transformation form is obtained as dmj
1
∂cave dm Ω = (Ω + 1) ∂θ dθ
The mass transfer term can be written as dmj
∫0
(17)
1 ∂ ⎛⎜ Ω ∂c ⎞⎟ j r r dr + r Ω ∂r ⎝ ∂r ⎠
1
(22)
Consequently the average concentration in the particle phase can be calculated by the Ωth moment as
Insert the right-hand side of eq 16, dmj
∂cave Ω r dr ∂θ
∂c dm Ω = ave ∂θ dθ
The basic idea of the moment method is to follow the time evolution of the moments of the concentration profiles, instead of modeling the concentration profiles directly. Then, the concentration profiles are calculated from the moment values by linear transformations.11 The remaining part of this section describes the general formulation of the moment method to solve the 1D model. The definition of the jth moment of ∂c/∂θ is the following (the subscript “i” representing individual component is left out of the equations for simplicity): dmj
1
Because ∂cave/∂θ is independent of the particle radius r, eq 22 is written as
∂c r = 0, i = 0; ∂r
∂ci = Bii(cbi − ci , r = 1) ∂r
∫0
(21)
This indicates that the Ωth moment describes the total mass change in the particle phase. In this study, the derivative of the first three moments (0th, 1st, and 2nd moments) are followed and the terms of ∂c/∂r, cc, and cw in eq 20 are obtained by the orthogonal polynomials that are functions of r2. The polynomial weights are calculated from these moments by a linear transformation.11
re = 4099
Yi , j − Yi , j ,3D Yi , j ,3D
× 100% (27)
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research
Article
3. RESULTS AND DISCUSSION 3.1. Model Parameters. The presented parameter values15,16 in Table 1 are intended to describe the natural Table 1. Wood Chip Parameters for the 3D Model dimensions (m) X Y Z
tortuosity τyz τxz τxy
0.004 0.01 0.02
porosity
2 3.7 3.7
εyz εxz εxy
0.7 0.7 0.7
characteristics of the wood chips, e.g. common dimensions and anisotropic properties of the wood chips, etc. In Table 2, it is Figure 1. Average concentration as a function of time in the particle phase. Dashed lines are the homogeneously mixed-chip phase (HMC) model. Solid lines are the 3D model with the finite difference method (3D_FD).
Table 2. General Parameters in the 3D, 1D, and Homogeneously Mixed-Chip Modelsa ICP (mol/L) C W P
ICB (mol/L)
0.36 0.566 0
Kr (l/mol·s) −3
−2.3 × 10 −2.3 × 10−3 +2.3 × 10−3
0.5 0 0
εb 0.65
a
ICP is the initial concentration in the particle phase; ICB is the initial concentration of the bulk phase.
shown that the wood chips are impregnated with “C” (representing, e.g., the alkali in the pulping process) in the beginning mimicking the fact that wood chips are typically fully impregnated with chemicals before pulping. The reaction rate and mass transfer coefficients17 are in Tables 2 and 3. According to Thiele modulus analysis,18 mass transfer resistance in this work is the rate-determining step. In this study, a number of different sets of parameters were tested and provided very similar results. Therefore only one set of parameters is listed in Tables 1−3. The simulation results are presented in the next section. The diffusion coefficient D and the external mass transfer coefficient k are estimated by using literature data.15 The values are shown in Table 3. For the homogeneously mixed-chip phase model, the diffusion coefficient D and external mass transfer coefficient k are lumped to the direction dependent parameters kxy, kxz, and kyz, which are calculated from the empirical correlations in eq 5. For the 1D model, the estimation methods of geometrical shape parameter Ω, the particle characteristic dimension L, and tortuosity τ were presented in section 2.3.1. 3.2. Simulation Results. The first pair of models to be compared is the 3D model (section 2.1) and the homogeneously mixed-chip phase model (section 2.2). The 3D model is solved by the finite difference method with nine grid points in each direction. In this work, the concentration profiles of 1/8 of a rectangular wood chip are solved. The complete profile in the whole wood chip is obtained by the symmetric feature of the rectangular chip. Therefore, the total number of grid points used to solve the model is 9 × 9 × 9. Figures 1 and 2 present
Figure 2. Average concentration as a function of time in the bulk phase. Dashed lines are the homogeneously mixed-chip phase model. Solid lines are the 3D model with the finite difference method (3D_FD).
the simulated average concentration profiles in particle phase and in bulk phase, respectively. It is obvious that there are major differences between the profiles produced by these two models. The quantitative deviation of the two models is 1.37 × 10−1, and the relative errors of most individual points are in the range of 3.0−40.0% as shown in Table 4. The highest relative error, about 50.0%, is for the component P in the particle phase. This indicates that ignoring the concentration profiles inside the chips may result in serious errors in the predicted compositions. The homogeneously mixed-chip phase model therefore is not applicable for accurate simulation although the number of variables and CPU time consumed are very low. For the 1D model (section 2.3), which is solved with the finite difference method, nine grid points are selected to guarantee the computational accuracy. The average concentrations in the particle phase and the bulk phase are presented in Figures 3 and 4. The difference between the 1D model and the 3D model becomes much smaller compared to the
Table 3. Parameters in the 3D and 1D, Homogeneously Mixed-Chip Phase, and 1D Models 3D and 1D C W P
HMC
D (m2/s)
k (m/s)
1.0 × 10−8 0 3.0 × 10−10
2.0 × 10−3 0 1.0 × 10−4
C W P
kxy (m/s)
kyz (m/s)
kxz (m/s)
4.05 × 10−6 0 1.21 × 10−7
7.47 × 10−6 0 2.24 × 10−7
4.05 × 10−6 0 1.21 × 10−7
4100
1D Ω L (m) τ
0.441 0.002 2
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research
Article
Table 4. Performance of Various Models and Numerical Schemesa model
numerical method
3D model HMC model 1D model 1D model
finite difference method finite difference method moment method
accuracy
variables
CPU time (s)
Dev
A-re
H-re
bad satisfatory satisfatory
9×9×9 1 9 3
96.1 0.04 2.4 0.3
1.37 × 10−1 7.87 × 10−3 8.29 × 10−3
3.0−40.0% 3.0% 3.0%
50.0% 5.6% 6.0%
a
Dev is the average quantitative deviation from the 3D model. A-re and H-re are the average relative error and the highest relative error of individual point. The computer applied in this study is with Intel Xeon CPU E31230 @ 3.20 GHz.
Figure 5. Average concentrations as a function of time in the particle phase. Dashed lines are the 1D model with the moment method (1D_MM). Solid lines are the 3D model with the finite difference method (3D_FD).
Figure 3. Average concentrations as a function of time in the particle phase. Dashed lines are the 1D model solved with the finite difference method (1D_FD). Solid lines are the 3D model with the finite difference method (3D_FD).
6. As can be seen, the profiles are also close to those obtained with the 3D model and are very similar to the profiles obtained
Figure 4. Average concentrations as a function of time in the bulk phase. Dashed lines are the 1D model solved with the finite difference method (1D_FD). Solid lines are the 3D model with the finite difference method (3D_FD).
Figure 6. Average concentrations as a function of time in the bulk phase. Dashed lines are the 1D model with the moment method (1D_MM). Solid lines are the 3D model with the finite difference method (3D_FD).
differences between the 3D and the homogeneously mixed chip phase models. The quantitative deviation between the two models is 7.87 × 10−3, and the relative errors of most individual points are less than 3.0% as reported in Table 4. The highest relative error, about 5.6% for component P, is in the particle phase. This indicates that the 1D model predictions are very close to the more complicated 3D model, which was not the case with the homogeneously mixed chip phase model. As shown in Table 4, the number of variables decreased dramatically compared to the 3D model. As a result, the CPU time needed to solve the 1D model decreased also significantly. The final comparison is between the two numerical methods to solve the 1D model, namely finite difference and moment method. The average concentration profiles of the 1D model solved with the moment method are presented in Figures 5 and
with the 1D model solved with the finite difference method (Because the profiles of 1D model solved with the moment method and the finite difference method overlap closely to each other, only the profiles with the moment method are presented in Figure 5 and 6). The quantitative deviation of the two models is 8.29 × 10−3, and the relative errors of most individual points are less than 3.0% as shown in Table 4. The highest relative error, about 6.0%, is for component P in the particle phase. As shown in Table 4, the number of variables and CPU time decreased further compared to the finite difference method. In this study, solving the 3D model involves 3 components with 9 × 9 × 9 grid points. This results in a system with 3 × 9 × 4101
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research
■
9 × 9 = 2187 coupled ordinary differential equations (ODEs). For the particle phase of the 1D model with nine grid points in the finite difference method, the number of ODEs reduces to 3 × 9 = 27. To solve the 1D model with the moment method, for each component, the 0th, 1st, and 2nd moments are conserved. This results in only 3 × 3 = 9 ODEs to be solved. Thus, the use of the 1D model and moment method significantly reduces the number of variables and saves greatly the computational power, which can be seen from the decreased CPU time reported in Table 4. This is especially beneficial when the system consists of tens of components and a great number of reactions, as well in cases where large scale digester units are modeled with appropriate material balances and fluid flow patterns. Furthermore, once the bulk phase model is employed, the number of ODEs in the 3D model and in the 1D model solved with finite difference method increases consequently. For the moment method, it is more convenient to calculate the average concentration of the bulk phase directly from the Ωth moment as in eq 25. Another attractive feature of the moment method is that the spatial PDE solution inherently conserves mass.11 In this work the moment method gives very low total mass balance error, which is close to the time integrator tolerance. As a comparison, the total mass balance error when using the finite difference method is, although within the acceptable limit, relatively high i.e. between 0.05 and 0.2%.
Article
AUTHOR INFORMATION
Corresponding Author
*E-mail: zheng.liu@aalto.fi. Notes
The authors declare no competing financial interest.
■ ■
ACKNOWLEDGMENTS Financial support from FIBIC Ltd (EffFibre and Fubio JR2 Research Programmes) is gratefully acknowledged.
4. CONCLUSIONS A general methodology to reduce the computational burden when modeling the mass transfer and reactions in anisotropic biomass particles was introduced. The pulping of wood chips was selected as a realistic application example to discuss the methodology. For modeling mass transfer and reactions in wood chips, three modeling approaches were compared in this article: 3D model (the most realistic, considered as the reference case), homogeneously mixed-chip phase model (traditionally used in digester models), and a novel 1D model. By comparing the results obtained with homogeneously mixed chip phase model (which ignores the concentration profiles inside the chips) to the results obtained with the 3D model, it was shown that performance of homogeneously mixed chip phase model was not adequate. However, the CPU time consumed to solve the 3D model was high. In this article, the equations for estimation of 1D model parameters, i.e., geometrical shape parameter, particle characteristic dimension, and a few other parameters were presented. The 3D model and the 1D model were solved with the finite difference method. The average concentration profiles inside the chips as a function of time obtained with the 1D model were reasonably close to those obtained with the 3D model. The number of variables to solve the 1D model was much less than in the 3D model. In order to reduce the number of variables and lessen the computational burden further, a high-order moment based weighted residual method was employed to solve the 1D model. The general formulation and implementation procedure of the moment method to solve the 1D model were introduced and discussed. Due to the greatly reduced number of ODEs, the moment method is capable of solving the 1D model relatively accurately with significantly decreased computational load, compared to the 3D model. This is especially beneficial when the system to be modeled consists of a great number of components and reactions.
NOMENCLATURE Aα = surface area of the particle [m2] Axy = surface area of coordinates x and y [m2] Axz = surface area of coordinate x and z [m2] Ayz = surface area of coordinate y and z [m2] Bi = Biot number, kiL/εpDpi [−] C = particle phase concentration [mol/m3] Cave = average concentration in the particle phase [mol/m3] Cb = liquid bulk phase concentration [mol/m3] Cref = concentration used for nondimensionalization [mol/ L3] c = dimensionless concentration [−] D = diffusion coefficient in water [m2/s] E′x = effective capillary cross-sectional area at x coordinate [−] Ey′ = effective capillary cross-sectional area at y coordinate [−] E′z = effective capillary cross-sectional area at z coordinate [−] Kr = reaction rate constant [m3/mol·s] k = external film diffusion coefficient [m/s] kr = dimensionless reaction kinetic coefficient, KrCrefL2/Dp_ref [−] kxy = lumped external film diffusion coefficient at the area of coordinate x and y, [m/s] kxz = lumped external film diffusion coefficient at the area of coordinate x and z [m s−1] kyz = lumped external film diffusion coefficient at the area of coordinate y and z [m/s] L = particle characteristic dimension [m] m = moment of a distribution or total number of simulation points [−] N = flux [mol/(s m2)] n = number of components [−] Sr = reaction term [mol/m3·s] R = coordinate of particle characteristic dimension [m] r = dimensionless coordinate of particle characteristic dimension [−] t = time [s] V = volume [m3] Vp = volume of the rectangular wood chip [m3] Vα = volume of the particle [m3] X = wood chip length of x coordinate [m] Y = wood chip length of y coordinate [m] Z = wood chip length of z coordinate [m]
Greek Symbols
Ω = geometrical shape parameter [−] θ = dimensionless time [−] η = dimensionless constant, Dpi/Dp_refτ [−] εb = bulk phase volume fraction in the bulk-particle phase system [−] εp = particle total porosity [−]
4102
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
Industrial & Engineering Chemistry Research
Article
εxy = porosity at the area of coordinate x and y [−] εxz = porosity at the area of coordinate x and z [−] εyz = porosity at the area of coordinate y and z [−] τ = particle tortuosity [−] τxy = tortuosity at the area of coordinate x and y [−] τxz = tortuosity at the area of coordinate x and z [−] τyz = tortuosity at the area of coordinate y and z [−]
(17) Sixta, H. Handbook of pulp: Chemical Pulp Process; Wiley: New York, 2008. (18) Hill, C. An Introduction to Chemical Engineering and Reactor Design; John Wiley & Sons, Inc.: Hoboken, NJ, 1977; pp 440−446.
Subscripts and Superscripts
ave = average b = liquid bulk phase i = component index j = jth moment or jth simulation point ref = reference value x = x coordinate or x direction of the wood chip y = y coordinate or y direction of the wood chip z = z coordinate or z direction of the wood chip
■
REFERENCES
(1) Ragauskas, A. J.; Williams, C. K.; Davison, B. H.; Britovsek, G.; Cairney, J.; Eckert, C. A.; Frederick, W. J., Jr.; Hallett, J. P.; Leak, D. J.; Liotta, C. L.; Mielenz, J. R.; Murphy, R.; Templer, R.; Tschaplinski, T. The path forward for biofuels and biomaterials. Science 2006, 311, 484−489. (2) Huang, H.-J.; Ramaswamy, S.; Tschirner, U. W.; Ramarao, B. V. A review of separation technologies in current and future biorefineries. Sep. Purif. Technol. 2008, 62, 1−21. (3) Gullichsen, J.; Fogelholm, C.-J. Papermaking science and Technology: 6. Chemical Pulping; Tappi Press: Finland, 2000. (4) Sjöström, E. Wood Chemistry: Fundamentals and Applications; Academic Press, 1993. (5) Simão, J. P. F.; Egas, A. P. V.; Carvalho, M. G.; Baptista, C. M. S. G.; Castro, J. A. A. M. Heterogeneous studies in pulping of wood: Modelling mass transfer of alkali. Chem. Eng. J. 2008, 139, 615−621. (6) Inalbon, M. C.; Mussati, M. C.; Mocchiutti, P.; Zanuttini, M. A. Modeling of Alkali impregnation of eucalyptus wood. Ind. Eng. Chem. Res. 2010, 50, 2898−2904. (7) Jutila, E. A. A. A Survey of Kraft Cooking Control Models. Instrumentation and Automation in the Paper, Rubber, Plastics and Polymerization Industries, Proceedings of the 4th IFAC Conference, Ghent, Belguim, June 3−5, 1980. (8) Grénman, H.; Wärnå, J.; Mikkola, J.-P.; Sifontes, V.; Fardim, P.; Murzin, D. Y.; Salmi, T. Modeling the influence of wood anisotropy and internal diffusion on delignification kinetics. Ind. Eng. Chem. Res. 2010, 49, 9703−9711. (9) Alopaeus, V.; Laavi, H.; Aittamaa, J. A dynamic model for plug flow reactor state profiles. Comput. Chem. Eng. 2008, 32, 1494−1506. (10) Roininen, J.; Alopaeus, V. The moment method for onedimensional dynamic reactor models with axial dispersion. Comput. Chem. Eng. 2011, 35 (3), 423−433. (11) Liu, Z.; Roininen, J.; Pulkkinen, I.; Sainio, T.; Alopaeus, V. Moment based weighted residual method − New numerical tool for a nonlinear multicomponent chromatographic general rate model. Comput. Chem. Eng. 2013, 53, 153−163. (12) Burghardt, A.; Kubaczka, A. Generalization of the effectiveness factor for any shape of a catalyst pellet. Chem. Eng. Process.: Process Intens. 1996, 35, 65−74. (13) Neretnieks, I. Analysis of some washing experiments of cooking chips. Svensk Papperstidning 1972, 75, 819−825. (14) Gu, T.; Tsai, G.-J.; Tsao, G. T. Modeling of nonlinear multicomponent chromatography. Adv. Biochem. Eng. Biotechnol. 1993, 49, 52−54. (15) Robertsen, L. Diffusion in wood: Part 4. Diffusion in radial and longitudinal direction. Ph.D. dissertation, Department of pulping technology, Åbo Akademi, Turku, Finland, 1993. (16) Stone, J. The Effective Capillary Cross-Sectional Area of Wood as a Function of pH. Tappi 1957, 40, 539. 4103
dx.doi.org/10.1021/ie403400n | Ind. Eng. Chem. Res. 2014, 53, 4096−4103
