Multi-Scale Modeling of the Dynamics of a Fibrous Reactor: Use of an Analytical Solution at the Micro-Scale to Avoid the Spatial Discretization of the Intra-Fiber Space

Abstract

Direct modeling of time-dependent transport and reactions in realistic heterogeneous systems, in a manner that considers the evolution of the quantities of interest in both, the macro-scale (suspending fluid) and the micro-scale (suspended particles), is currently well beyond the capabilities of modern supercomputing. This is understandable, since even a simple system such as this can easily contain over 10⁷ particles, whose length and time scales differ from those of the macro-scale by several orders of magnitude. While much can be gained by applying direct numerical solution to representative model systems, the direct approach is impractical when the performance of large, realistic systems is to be modeled. In this study we derive and analyze a “hybrid” model that is suitable for fibrous reactors. The model considers convection/diffusion in the bulk liquid, as well as intra-fiber diffusion and reaction. The essence of our approach is that diffusion and (first-order) reaction in the intra-fiber space are handled semi-analytically, based on well-established theory. As a result, the problem of intra-fiber transport and reaction is reduced to an easily solvable set of $n_0$ ODEs, where $n_0$ is the number of terms in the Bessel expansion evaluated without recourse to approximation; this set is coupled, point-wise, with a numerical model of the macro-scale. When the latter is discretized using N nodes, the total “hybrid” model for the system consists of a system of $N\left(2+n_0\right)$ ODEs, which is easily solvable on a modest workstation. Parametric analyses are presented and discussed.

1. Introduction

The modeling of coupled flow and transport/reaction in systems comprised of a solid phase dispersed in a fluid phase has long been of interest to scientists and engineers. These systems involve heat/mass transfer processes occurring at multiple time/length scales [1,2,3] and are found in a wide range practical applications. Examples range from bioreactors [4] to catalytic cracking [5], where the solid phase may be mobile, as in fluidized bed reactors, or fixed in place. Other areas of application are found in the food industry, in the form of drying [6] and meat processing [7]. The particles themselves can be porous, so a complete model must take into account transport in the carrier fluid, as well as the exchange with and diffusion/reaction within the particles.

A single complete analytical solution is not possible for such systems, so many different approaches have been used in their simplification and solution. The problem can be simplified by focusing on the steady state solution and, in some cases, ignoring the diffusion within the particles [8]. However, even these assumptions still lead to a demanding numerical problem. In addition, transient solutions are desirable as they allow more information to be extracted, especially when model predictions are compared with dynamic experiments [6]. Another approach is to focus on only a single spatial scale, modeling diffusion in a single structure [7]. Internal diffusion has also been taken into account in countercurrent plug flow problems by relating the average concentration within the particles to the local concentration in the fluid [9].

In the past few decades, advances in multi-scale and homogenization theories [10], as well as in Computational Fluid Dynamics (CFD) codes and computing facilities, have pushed the boundaries of what can be solved numerically [11]. Particle Resolved Simulations can solve the local heat/mass transfer equations at the particle scale for both fixed and fluidized beds [12]. Detailed particle-based approaches have been proposed using the simulation of a small number of particles in repeating unit cells [3]. One drawback of directly solving the microscale transport equations in every single particle is that prohibitively large meshes are required for, macroscopically, large systems. The problem is more acute in time-dependent simulations, and poses a severe restriction in simulating real-life situations, in which micron- or sub-micron-sized particles may be used.

The purpose of this work is to derive and present a hybrid (semi-analytical) solution to the problem of coupled micro/macro scale mass-transfer in reactors containing fibers (micro-scale) oriented transversely to the bulk flow. Besides intra-fiber diffusion, we also consider a first-order reaction to be taking place within the fibers. This work utilizes the approach previously used for reactors containing spherical catalyst particles and extends it to fibrous systems [13]. The key assumption in the solution derivation is that individual fibers are treated as if they are surrounded by a liquid of a uniform albeit time-varying concentration. Besides a small fiber diameter, this assumption is more likely to hold when the bulk flow is transverse to the fiber axis, implying that our systems also consist of unidirectional fibers. The concentration in the liquid phase (macro-scale) is resolved using the well-established dispersion model. The proposed model only requires meshing of the macro-scale, while the intra-particle concentration is described by a semi-analytical formula which requires the solution of a small number of ODEs. This is easily tractable using standard ODE solvers.

Although this specific study is using the axial dispersion model to describe the transport in the bulk liquid, the fiber-scale model is portable to any other macroscopic situation, as long as the fibers are small relative to the concentration gradients in the bulk liquid. In Section 2, the development of the model will be outlined for the case of constant reactor porosity, fiber porosity, and radius; this is not a necessary assumption and the model can be modified to accept spatially varying parameters.

2. Model Development

The general approach of the outlined semi-analytical solution is to use the local surface concentration ($C_R$) in the fiber phase as a time-varying boundary condition for the intra-fiber transport problem. The analytical solution for the concentration profile in the fiber is approximated by expanding the solution in Bessel functions. The concentration profile allows for the calculation of the flux from the liquid phase to the fiber phase, yielding an equation for the local liquid concentration ($C_L$). Discretization of the bulk scale leads to a system of ODEs which can be easily solved using standard numerical methods.

2.1. Problem Description

In order to describe the dynamic response of the system shown in Figure 1, the governing differential equations for the mass transfer in the bulk liquid and the mass transfer/consumption in the intra-fiber space must be solved. The boundary condition at the surface of the fiber governs the mass transfer between the fibers and the bulk liquid, therefore providing the link between the micro (intrafiber spaces) and macro (bulk liquid) scales of the problem. The axial distance along the reactor and the radial co-ordinate in a fiber are the two spatial independent variables, while the time is the third independent variable. In a manner similar to earlier studies [13,14] the following phenomena are taken into account in the following analysis:

• Mass accumulation in the bulk liquid.
• Mass diffusion from the bulk liquid, across the boundary layer, and into the fibers.
• Mass diffusion, accumulation and consumption (via a first-order chemical reaction) in the intra-fiber space.
• Axial dispersion in the bulk liquid.

The operating assumptions of the model are:

• The fiber locations are fixed.
• The system is isothermal.
• The first order reaction constant and the intra-fiber diffusion coefficient are constant.
• The bulk liquid concentration gradients in the perpendicular plane (along the fiber length), fiber aging and the pressure drop in the reactor are negligible. This implies that the bulk flow is strictly transverse to the fiber axes.
• The fibers are of uniform size and uniformly dispersed along the reactor.

2.2. Model Development

The following equation describes the mass transfer, accumulation, and consumption in the intra-fiber space:

$${\epsilon }_f\frac{\partial C_f}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{D}_{eff}\frac{\partial C_f}{\partial r}\right)-{\epsilon }_{f}k{C}_f$$

(1)

where $C_f$ is the concentration in the intra-fiber pore volume, $D_{eff}$ is the effective intrafiber diffusivity, $k$ is the first order reaction rate constant, and ${\epsilon }_f$ is the fiber void fraction. The initial and boundary conditions are given by

$$C_f\left(r,t=0\right)=0$$

(2)

$${\left(\frac{\partial C_f}{\partial r}\right)}_{r=0}=0$$

(3)

$$-D_{eff}{\left(\frac{\partial C_f}{\partial r}\right)}_{r=R}=k_e\left(C_R^L-C_L\right)$$

(4)

where $C_L$ is the concentration in the bulk liquid, $C_R^L$ is the liquid concentration at the inner side of the boundary layer around the fiber, and $k_e$ is the external mass transfer coefficient. $C_R^L$ is related to the fiber surface concentration ($C_R$) via a partition coefficient, $\gamma =C_R^L/C_R$. The bulk-fluid concentration C_L Equation (4) links the micro- and the macro-scale processes. This microscale flux across the fiber surface Equation (4) can be multiplied by the total fiber surface area per unit volume of reactor to yield the flux from the bulk liquid to the fibers at any axial position in the reactor, $Q^{*}$.

$$Q^{*}=2\left(1-\epsilon \right)\left(\frac{D_{eff}}{R\left(F/V_T\right)}\right){\left(\frac{\partial C_f}{\partial r}\right)}_{r=R}$$

(5)

where $\epsilon$ is the reactor void fraction, $R$ is the fiber radius, $F$ is the bulk flow rate, and $V_T$ the total reactor volume.

Using the dispersion model, in dimensionless form, the mass balance in the bulk liquid is

$$\frac{d{C}_L}{d{t}_R}=\frac{1}{Pe}\frac{{\partial }^2{C}_L}{\partial x^2}-\frac{\partial C_L}{\partial x}-Q^{*}$$

(6)

in which, $Pe$ is the Péclet number, describing the relative strength of convection and diffusion in the bulk liquid, $x$ is the dimensionless axial position, and $t_R$ is the dimensionless time. Evidently, the flux term ($Q^{*}$) in Equation (5) links the micro Equation (1) and macro Equation (6) scales of the problem. The three terms on the right-hand side of Equation (6) are the mass flux due to diffusion, convection, and transfer to the fibers. The initial and boundary conditions are

$$C_L\left(x,0\right)=0$$

(7)

$${\left(\frac{\partial C_L}{\partial x}\right)}_{x=0}=-Pe\left(C_0-C_L\right)$$

(8)

$${\left(\frac{\partial C_L}{\partial x}\right)}_{x=1}=0$$

(9)

As explained in Appendix A, the analytical solution of Equation (1) (available at constant surface concentration) is used to derive an expression for the intra-fiber concentration, subject to a time-varying surface concentration. The result is that the fiber surface concentration ($C_R$) can be determined by solving the following ODE:

$$\frac{\partial C_R}{\partial t_R}=\frac{B_m\epsilon D_R}{2{\epsilon }_f{S}_i}{C}_L-\left(\frac{B_m\epsilon D_R}{2{\epsilon }_f{S}_i}\gamma +9{\varphi }^2\epsilon D_R\right)C_R-\frac{\epsilon D_R}{S_i}{\displaystyle \sum }_{n=1}^{n_0}{\mathsf{\Psi }}_{\mathrm{n}}$$

(10)

The ${\mathsf{\Psi }}_{\mathrm{n}}$ terms appearing in Equation (11) are described by ODEs of the form.

$$\frac{\partial {\mathsf{\Psi }}_{\mathrm{n}}}{\partial t_R}+\epsilon D_R\left[{\beta }_n{}^2+9{\mathsf{\Phi }}^2\right]{\mathsf{\Psi }}_{\mathrm{n}}=\frac{\partial C_R}{\partial t_R}+9\epsilon {\mathsf{\Phi }}^2{D}_R{C}_R$$

(11)

As a result, and as explained in the Appendix A, the fluid-to-fiber mass flux can be evaluated either via the difference between the surface and liquid concentration Equation (12) or using the summation of the ${\mathsf{\Psi }}_{\mathrm{n}}$ terms Equation (13).

$$Q^{*}=-2\left(1-\epsilon \right)D_R{B}_m\left(\gamma C_R-C_L\right)$$

(12)

$$Q^{*}=4\left(1-\epsilon \right)D_R{\displaystyle \sum }_{n=1}^{\infty }{\mathsf{\Psi }}_{\mathrm{n}}$$

(13)

In short, the microscale discretization of the fibers has been avoided and, instead, a system of $\left(2+n_0\right)$ ODEs describe the evolution of $C_L$, $C_R$, at a given axial position. Equation (6) can then be discretized in the axial direction and written, in dimensionless form at axial node $i$, as:

$$\frac{d{C}_{L,i}}{d{t}_R}=\frac{1}{Pe}\frac{C_{L,i+1}-2C_{L,i}+C_{L,i-1}}{\mathsf{\Delta }{x}^2}-\frac{C_{L,i+1}-C_{L,i-1}}{2\mathsf{\Delta }x}-2D_R{B}_m\left(1-\epsilon \right)\left(C_{L,i}-\gamma C_{R,i}\right)$$

(14)

The ODEs defined by Equations (10), (11), and (14) are applicable at each node, leading to a system of $N\left(2+n_0\right)$ ODEs, which can be solved using standard ODE solvers such as ode45 or ode113 available in MATLAB.

3. Results

We have performed a series of simulations using MATLAB ODE solvers ode45 or ode113 over a range of dimensionless parameter values. The Thiele Modulus, $\Phi$, describing the ratio of the reaction rate to the diffusion rate inside the fiber, varied from 0 to 100. The Biot number, $B_m$, describes the ratio of the external to internal mass transfer resistances varied from 10⁻³ to 100. The dimensionless diffusivity, $D_R$, was varied from 0.001 to 100. Fiber voidage, ${\epsilon }_f$ and Péclet number, $Pe$, were also varied.

Model predictions concerning the effect of the Péclet number, in the presence and absence of intra-fiber porosity, are shown in Figure 2. When porous fibers are considered, there is transfer of the reactant from the bulk as the species diffuse into the fibers. This results in a modified effluent profilecompared to the impermeable (solid) fibers case. This effect is more pronounced at high Péclet numbers.

The effect of the Biot number is shown in Figure 3 for two different levels of diffusivity for a non-reacting species ($\Phi$ = 0). In each case, very low Biot numbers lead to the condition where the liquid concentration profile is limited by a high external mass transfer resistance to the fibers, leading to a negligible amount of exchange between the two phases. On the other end of the spectrum the internal resistance of the fibers dominates. The Biot numbers that correspond to these limiting cases are dependent on the other system parameters, such as $D_R$. For the intermediary cases, the $B_m$ and $D_R$ terms both play a role as neither internal, nor external mass transfer resistances dominate.

In Figure 4, the Thiele modulus is varied from 0 to 100 and the effect on the effluent profile and the radial concentration profile, computed from Equation (A5) within a fiber at the end of the reactor at $t_R=10$ are shown. When $\Phi$ is very low, the reaction is rate-limited as the reactant is consumed very slowly and builds up in the fiber. When $\Phi$ is very large the reaction in the fibers proceeds very quickly and the mass flux to the fiber is not large enough to sustain a high concentration in the fiber; this corresponds to the mass transfer limited regime. In between the two extremes, all the parameters play their role in the process.

4. Case Study

A hybrid, semi-analytical model for coupled flow and mass-transfer processes in reactors containing fibers oriented perpendicular to the direction of bulk flow has been presented in Section 2, with some results of parametric analyses being shown in Section 3. This model is also applicable to situations where the fiber radii/porosities vary along the reactor. This is because the relationship we derive between the surface concentration in the fiber and that of the bulk liquid apply at each axial position. The appropriate constants from Equations (11)–(13) would simply need to be replaced by their values at node $i$. Another example would be the case where the reactor voidage is not constant along the reactor, representing bundles of fibers placed at different locations. Such a case is shown in Figure 5, in which a decreased voidage over 10 percent of the reactor length is introduced at three positions: near the inlet, at the center, and near the outlet. The constant voidage example having the same average $\epsilon$ over the reactor is also shown.

The results of the analysis of the cases represented in Figure 5, in terms of effluent concentrations and axial concentration profiles at steady state, are shown in Figure 6. With all other factors remaining constant, the location of the fiber bundle does appear to slightly affect the transient response of the system (Figure 6a), the fiber bundles leading to a slightly increased conversion of the feed material within the reactor. This example illustrates one of the advantages of the semi-analytical model, as it can easily capture differences that would have been ignored by the simplifications used in a model based on bulk/effective properties.