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

: 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 ﬂuid) 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 7 particles, whose length and time scales di ﬀ er 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 ﬁbrous reactors. The model considers convection / di ﬀ usion in the bulk liquid, as well as intra-ﬁber di ﬀ usion and reaction. The essence of our approach is that di ﬀ usion and (ﬁrst-order) reaction in the intra-ﬁber space are handled semi-analytically, based on well-established theory. As a result, the problem of intra-ﬁber 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 ( 2 + n 0 ) ODEs, which is easily solvable on a modest workstation. Parametric analyses are presented and discussed. for a dependent


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,

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.

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 Fluids 2020, 5, 3 3 of 11 third independent variable. In a manner similar to earlier studies [13,14] the following phenomena are taken into account in the following analysis: i.
Mass accumulation in the bulk liquid. ii.
Mass diffusion from the bulk liquid, across the boundary layer, and into the fibers. iii.
Mass diffusion, accumulation and consumption (via a first-order chemical reaction) in the intra-fiber space. iv.
Axial dispersion in the bulk liquid.
i. Mass accumulation in the bulk liquid. ii. Mass diffusion from the bulk liquid, across the boundary layer, and into the fibers. iii. Mass diffusion, accumulation and consumption (via a first-order chemical reaction) in the intra-fiber space. iv. Axial dispersion in the bulk liquid.
The operating assumptions of the model are: i.
The fiber locations are fixed. ii. The system is isothermal. iii. The first order reaction constant and the intra-fiber diffusion coefficient are constant. iv. 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. v. The fibers are of uniform size and uniformly dispersed along the reactor.

Model Development
The following equation describes the mass transfer, accumulation, and consumption in the intra-fiber space: where is the concentration in the intra-fiber pore volume, is the effective intrafiber diffusivity, is the first order reaction rate constant, and is the fiber void fraction. The initial and boundary conditions are given by where is the concentration in the bulk liquid, is the liquid concentration at the inner side of the boundary layer around the fiber, and is the external mass transfer coefficient. is related to the fiber surface concentration ( ) via a partition coefficient, = / . The bulk-fluid concentration CL 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, * . * = 2(1 − ) ( / ) (5) The operating assumptions of the model are: i. The fiber locations are fixed. ii.
The system is isothermal. iii.
The first order reaction constant and the intra-fiber diffusion coefficient are constant. iv.
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. v.
The fibers are of uniform size and uniformly dispersed along the reactor.

Model Development
The following equation describes the mass transfer, accumulation, and consumption in the intra-fiber space: where C f is the concentration in the intra-fiber pore volume, D e f f is the effective intrafiber diffusivity, k is the first order reaction rate constant, and ε f is the fiber void fraction. The initial and boundary conditions are given by − D e f f ∂C f where C L is the concentration in the bulk liquid, C L R 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 L R is related to the fiber surface concentration (C R ) via a partition coefficient, γ = C L R /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 * . where ε 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 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 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: The Ψ n terms appearing in Equation (11) are described by ODEs of the form.
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 Ψ n terms Equation (13).
In short, the microscale discretization of the fibers has been avoided and, instead, a system of (2 + n 0 ) 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: The ODEs defined by Equations (10), (11), and (14) are applicable at each node, leading to a system of N(2 + n 0 ) ODEs, which can be solved using standard ODE solvers such as ode45 or ode113 available in MATLAB.

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, Φ, 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 −3 to 100. The dimensionless diffusivity, D R , was varied from 0.001 to 100. Fiber voidage, ε 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.
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 ( = 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 . For the intermediary cases, the and terms both play a role as neither internal, nor external mass transfer resistances dominate.
The effect of the Biot number is shown in Figure 3 for two different levels of diffusivity for a non-reacting species (Φ = 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.
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 ( = 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 . For the intermediary cases, the and terms both play a role as neither internal, nor external mass transfer resistances dominate.
(a) (b)   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 Φ is very low, the reaction is rate-limited as the reactant is consumed very slowly and builds up in the fiber. When Φ 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. 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 = 10 are shown. When is very low, the reaction is rate-limited as the reactant is consumed very slowly and builds up in the fiber. When 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.

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 . 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 over the reactor is also shown.

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 ε over the reactor is also shown. 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 = 10 are shown. When is very low, the reaction is rate-limited as the reactant is consumed very slowly and builds up in the fiber. When 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.

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 . 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 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. 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.  Funding: This work was carried out with support from Grant Award Number 090118FD5313 under the "Structure-Property Correlations in Multi-Scale Composites" project at Nazarbayev University.

Conflicts of Interest:
The authors declare no conflict of interest.

Notation of Parameters
β n positive roots of J 0 (β) = 0, where J ν is a Bessel function of order ν ε reactor void fraction ε f intra-fiber void fraction γ partition coefficient, C L /C f Ψ n functions defined by Equations (6)

Intra-Fiber Concentration Profile in the Presence of Diffusion and 1st Order Reaction under a Time-Varying Surface Concentration
The first step in the solution to Equation (1) subject to a time-varying surface concentration relies on existing analytical solution for diffusion inside a cylinder [8,9] initially at zero concentration and subject to a constant surface concentration of C 0 for t > 0. This solution is: J ν is the Bessel function of order ν, β n are the positive roots of J 0 (β) = 0 and k n = D e f f β n R

2
. When a first-order reaction is present Danckwert's method [15] can be used to show that the time-dependent concentration profile in the fiber (C 2 ) for the same constant surface concentration is given by Using Equation (A1) to evaluate Equation (A2) leads to the following intra-fiber concentration profile for constant surface concentration and in the presence of diffusion and first-order reaction The normalized concentration profile under a time-varying surface concentration (C R ) can now be found using Duhamel's Theorem [8] Carrying out the integrations, the resulting expression for the intra-fiber profile, under a time-varying surface concentration, C R , is given by where the Ψ n functions are given by Fluids 2020, 5, 3 9 of 11 or in differential form by with initial conditions Ψ n (0) = 0. The derivative in Equation (5) can now be evaluated, using Equation (A5), as The relationship between the surface concentration and liquid concentration can therefore be written as And the flux term can be calculated either via the difference between the surface and the liquid concentration Equation (A10) or via the summation of the Ψ n terms Equation (A11).
The large/increasing values of the factor multiplying the Ψ n term (β n are monotonously increasing) in Equation (A7) imply that, at large n, the Ψ n functions follow the forcing function on the right hand side of the equation. Above a certain threshold, the Ψ n terms can be approximated using the quasi-steady state approximation The infinite summation in Equation (A9) is therefore split into two parts, with the first n 0 terms being evaluated using Equation (A7) and the rest of the summation to infinity being approximated by Equation (A12).
Examples of the time-varying Ψ n terms of different orders and their steady state approximations using Equation (A12) (dashed lines, denoted "S.S.") are shown in Figure A1a,b. The lower order n terms ( Figure A1a) are poorly approximated by Equation (A12) and clearly outweigh the higher order terms, which can be reasonably approximated by Equation (A12) ( Figure A1b). Computing more than the first ≈ 50 Ψ n terms does not noticeably impact the result of the summation for the example under consideration in Figure A1c. In general, more Ψ n terms will be necessary for a good approximation when the ∂C R ∂t R term is large. Computing more than the first ≈ 50 terms does not noticeably impact the result of the summation for the example under consideration in Figure A1c. In general, more terms will be necessary for a good approximation when the term is large. The relationship between the accuracy and is more clearly shown in Figure A2a where the maximum error (over the timescale shown in Figure A1) relative to the case of = 200 is shown for two different variables; namely, , and , . The number of Ψ terms to be calculated to reach a desired level of accuracy for a given reactor can further depend on the variable of interest, for example the effluent concentration requires fewer terms than the concentration at the fiber center to reach a certain accuracy threshold. As increases, the computation time also increases significantly, as shown in Figure A2b, therefore the should be selected to balance accuracy and computation time for each reactor. A model based on a similar method of solution of the multi-scale problem has been used in [16] for the prediction of the response of a fixed bed made of Ca-Alginate gel particles. It was found that model predictions were in excellent agreement to experimental results, which exhibited the same, qualitatively distinct, features predicted by the model.  The relationship between the accuracy and n 0 is more clearly shown in Figure A2a where the maximum error (over the timescale shown in Figure A1) relative to the case of n 0 = 200 is shown for two different variables; namely, C L, x = 1 and C f , r = 0 and x = 0 . The number of Ψ n terms to be calculated to reach a desired level of accuracy for a given reactor can further depend on the variable of interest, for example the effluent concentration requires fewer Ψ n terms than the concentration at the fiber center to reach a certain accuracy threshold. As n 0 increases, the computation time also increases significantly, as shown in Figure A2b, therefore the n 0 should be selected to balance accuracy and computation time for each reactor. A model based on a similar method of solution of the multi-scale problem has been used in [16] for the prediction of the response of a fixed bed made of Ca-Alginate gel particles. It was found that model predictions were in excellent agreement to experimental results, which exhibited the same, qualitatively distinct, features predicted by the model.
for the example under consideration in Figure A1c. In general, more terms will be necessary for a good approximation when the term is large. The relationship between the accuracy and is more clearly shown in Figure A2a where the maximum error (over the timescale shown in Figure A1) relative to the case of = 200 is shown for two different variables; namely, , and , . The number of Ψ terms to be calculated to reach a desired level of accuracy for a given reactor can further depend on the variable of interest, for example the effluent concentration requires fewer terms than the concentration at the fiber center to reach a certain accuracy threshold. As increases, the computation time also increases significantly, as shown in Figure A2b, therefore the should be selected to balance accuracy and computation time for each reactor. A model based on a similar method of solution of the multi-scale problem has been used in [16] for the prediction of the response of a fixed bed made of Ca-Alginate gel particles. It was found that model predictions were in excellent agreement to experimental results, which exhibited the same, qualitatively distinct, features predicted by the model.  Equations (A9) and (A13) combine to yield the differential equation for the surface concentration as a function of the bulk liquid concentration and the process parameters. where S i is the summation