Simulation of Convection – Diffusion Transport in a Laminar Flow Past a Row of Parallel Absorbing Fibers

A numerical simulation of the laminar flow field and convection–diffusion mass transfer in a regular system of parallel fully absorbing fibers for the range of Reynolds numbers up to Re = 300 is performed. An isolated row of equidistant circular fibers arranged normally to the external flow is considered as the simplest model for a hollow-fiber membrane contactor. The drag forces acting on the fibers with dependence on Re and on the ratio of the fiber diameter to the distance between the fiber axes, as well as the fiber Sherwood number versus Re and the Schmidt number, Sc, are calculated. A nonlinear regression formula is proposed for calculating the fiber drag force versus Re in a wide range of the interfiber distances. It is shown that the Natanson formula for the fiber Sherwood number as a function of the fiber drag force, Re, and Sc, which was originally derived in the limit of high Peclet numbers, is applicable for small and intermediate Reynolds numbers; intermediate and large Peclet numbers, where Pe = Re × Sc; and for sparse and moderately dense rows of fibers.


Introduction
Studies of convection-diffusion transfer in fibrous media are essential for the development of various technological processes.The results obtained are employed in solutions to the problems concerning the filtration of suspended particles [1,2], sorption, catalysis, electrochemistry, and especially, in membrane separation technology [3,4], where membranes in the form of hollow fibers are used for membrane gas absorption/desorption [5,6], ultra- [7] and nanofiltration [8], and gas separation [9,10].Within the field of membrane technology, a great number of studies is devoted to the modelling of different mass-transfer processes in hollow-fiber membrane contactors [11][12][13][14][15][16].When simulating the membrane separation of mixtures in gas-liquid hollow-fiber contactors, it is important to calculate the arrival of the solute at the streamlined fiber surface; that is, to determine the collected fraction of the solute (Sherwood number), which depends on the solute diffusion coefficient, the fiber diameter, and the flow conditions-the flow velocity, viscosity, and temperature.It should be noted that the solute transport towards the fiber surface is strongly influenced by the velocity distribution near the fiber.At low flow velocities, corresponding to low Reynolds numbers Re < 1, the flow field depends on the single parameter α-the fiber packing density, whereas at higher velocities, at Re ≥ 1, it is governed by the two parameters of α and Re.In studying transport processes in real porous media, the models with known flow fields are used in order to eliminate the influence of structural uncertainty.The ordered systems of parallel cylinders with circular cross sections, arranged normally to the laminar flow with square or hexagonal packing or a single row of parallel fibers are usually used as models [1].The so-called cell models with axially symmetric flow past a fiber in a cell are commonly used at low Reynolds numbers.However, they are nonapplicable for the case of Re > 1.They also were shown to be not well-suited for dense systems of fibers, such as hollow-fiber membrane contactors, which are characterized by a high packing density with α > 0.3 (a/h > 0.5, where a is the fiber radius and 2h is the distance between the axes of neighboring fibers).However, using the cell model, the role of the mutual hydrodynamic influence of fibers on the hydrodynamic resistance of the fiber system and on the diffusion mass transfer was clarified [1].
The mass and heat transfer have been studied in detail [1,[16][17][18][19][20][21][22][23] on the basis of analytical and numerical solutions for the flow fields obtained for the fiber lattices [24][25][26][27][28].As it turned out, a single row of fibers is a particularly convenient model when studying fluid flow and mass transfer at Re > 1 [16,20], when the flow symmetry is broken with increasing Re, and when vortices are formed behind the fibers, leading to an increase in the diffusion flux to the fibers.At Re > 1, the fiber drag force and the Sherwood number should increase.This increase is more noticeable when the system is more porous.Calculations for the fiber drag force and the fiber Sherwood number in a single row were performed in our previous work [16] for Re < 50 and for a single prescribed value of the Schmidt number, Sc.For this example, the fiber Sherwood number was found as a function of the inlet flow velocity U.The case considered in [16] corresponds to typical operating conditions of a membrane contactor employed for the separation of gases and liquids.However, in a number of processes, for example in membrane distillation, the Reynolds numbers can be several times higher.In this paper, we consider the diffusion mass transfer in a wider range of Reynolds numbers up to Re = 300.Based on the obtained flow field, the fiber drag force and the fiber retention efficiency are obtained versus the process conditions and the interfiber distance.

Flow Field in a System of Parallel Fibers at Re > 1
Let us consider the viscous incompressible flow at low and intermediate Reynolds numbers in a single row of parallel equidistant fibers arranged normally to the flow direction.The corresponding simulation cell is shown in Figure 1.density, whereas at higher velocities, at Re  1, it is governed by the two parameters of  and Re.In studying transport processes in real porous media, the models with known flow fields are used in order to eliminate the influence of structural uncertainty.The ordered systems of parallel cylinders with circular cross sections, arranged normally to the laminar flow with square or hexagonal packing or a single row of parallel fibers are usually used as models [1].The so-called cell models with axially symmetric flow past a fiber in a cell are commonly used at low Reynolds numbers.However, they are nonapplicable for the case of Re > 1.They also were shown to be not well-suited for dense systems of fibers, such as hollow-fiber membrane contactors, which are characterized by a high packing density with  > 0.3 ( h a > 0.5, where a is the fiber radius and h 2 is the distance between the axes of neighboring fibers).However, using the cell model, the role of the mutual hydrodynamic influence of fibers on the hydrodynamic resistance of the fiber system and on the diffusion mass transfer was clarified [1].
The mass and heat transfer have been studied in detail [1,1623] on the basis of analytical and numerical solutions for the flow fields obtained for the fiber lattices [2428].As it turned out, a single row of fibers is a particularly convenient model when studying fluid flow and mass transfer at Re > 1 [16,20], when the flow symmetry is broken with increasing Re, and when vortices are formed behind the fibers, leading to an increase in the diffusion flux to the fibers.At Re > 1, the fiber drag force and the Sherwood number should increase.This increase is more noticeable when the system is more porous.Calculations for the fiber drag force and the fiber Sherwood number in a single row were performed in our previous work [16] for Re < 50 and for a single prescribed value of the Schmidt number, Sc.For this example, the fiber Sherwood number was found as a function of the inlet flow velocity U .The case considered in [16] corresponds to typical operating conditions of a membrane contactor employed for the separation of gases and liquids.However, in a number of processes, for example in membrane distillation, the Reynolds numbers can be several times higher.In this paper, we consider the diffusion mass transfer in a wider range of Reynolds numbers up to Re = 300.Based on the obtained flow field, the fiber drag force and the fiber retention efficiency are obtained versus the process conditions and the interfiber distance.

Flow Field in a System of Parallel Fibers at Re > 1
Let us consider the viscous incompressible flow at low and intermediate Reynolds numbers in a single row of parallel equidistant fibers arranged normally to the flow direction.The corresponding simulation cell is shown in Figure 1.To find the flow field, we will numerically solve the Navier-Stokes equations in the steady-state approximation: where is the Reynolds number based on the fiber diameter,  is the fluid density,  is the dynamic viscosity,  is the Laplace operator,  is the nabla operator, is the flow To find the flow field, we will numerically solve the Navier-Stokes equations in the steady-state approximation: where Re = 2aUρ/µ is the Reynolds number based on the fiber diameter, ρ is the fluid density, µ is the dynamic viscosity, ∆ is the Laplace operator, ∇ is the nabla operator, u = {u, v} is the flow velocity vector, p = p * a/Uµ is the pressure, the symbol • is the scalar product, and * denotes the dimensional variable.Hereinafter, all the values are reduced to dimensionless form by normalization to the fiber radius a and to the inlet unperturbed velocity U.The no-slip condition u = 0 is set at the fiber surface 4, while the condition of undisturbed flow u = {1,0} is set at the entrance of the simulation cell 1; the condition of vanishing viscous stresses (zero gradients) is applied at the outlet boundary 3; and at the boundaries labelled 2, the conditions of symmetry for the velocity components are used.The dimensionless fiber drag force per unit length was found as the surface integral of the projection of the local total stress on the flow direction: where T = (−pI + σ )n is the local total stress, σ is the viscous stress tensor, I is the unit tensor, n is the outer normal vector to the surface, dΣ is the surface element, and S g is the fiber surface area.
The fiber drag force is related with the dimensionless pressure drop across the fiber row as shown by ∆p = Fa/2h.Figures 2-4 show the results of simulations of the fiber drag force in a row of fibers dependent on Re.It is seen from Figures 3 and 4 that at the limit of small Reynolds numbers, the curves tend asymptotically to straight lines, corresponding to the Stokes flow regime at small but nonzero Re numbers.The Stokes flow region at Re 1 is described at a/h < 0.5 by the analytical formula of Miyagi [24]: and at a/h > 0.7, it is governed by the formula derived in the lubrication approximation by Keller [29]: Fibers 2018, 6, x FOR PEER REVIEW 5 of 11 It should be noted that the F (Re) curve for an isolated cylinder does not reach the plateau at Re << 1.At the same time, the fiber drag force in the fiber row increases more slowly with Re, as compared with that of an isolated cylinder.This is explained by the mutual hydrodynamic influence of the fibers.In the next Section, the calculated flow field is used for modeling the diffusion transport in the row of fibers.

Convection-Diffusion Mass Transfer in a System of Absorbing Fibers
To determine the distribution of the solute concentration С in the flow, the elliptic equation of convection-diffusion is solved: where Schmidt number, D is the solute diffusion coefficient, and . Here, all the values are reduced to dimensionless form by normalizing to a , U , and the inlet concentration , is set at the fiber surface 4. At the entrance of the cell 1, the condition of uniform concentration 1  C is set, and at the outlet boundary 3, the zero-gradient condition, . At the boundaries 2, the conditions of symmetry are used.The method of solving The regression formula for estimating F in the whole range of a/h was proposed in [16].Our results of computations for the drag force on a fiber in a row and in lattices of fibers are in excellent agreement with the known experimental data, with the results of computations of other authors, and with other analytical formulas [1,[24][25][26]30].
Figures 3 and 4 show the onset of the flow inertia effect.It is seen that the denser the row, the later the onset of nonlinear inertial effects; that is, the higher the Re number is at which the drag force begins to increase.This effect was shown for the square array of cylinders in [31], where it was noted that the relative importance of inertia becomes smaller as the volume fraction approaches close packing, because the largest contribution to the dissipation in this limit comes from the viscous lubrication flow in the small gaps between the cylinders.Furthermore, in the region 100 ≤ Re ≤ 300, these curves, as can be seen from Figure 2, are linear and parallel, which allows us to fit them by a simple formula (with the mean average percentage error <ε> less than 0.6%): where Fibers 2018, 6, x FOR PEER REVIEW 4 of 11 also plotted (in Figures 2 and 3) the curves for an isolated cylinder by the piecewise continuous approximation formula from [32]: computed with the step 0.05 starting from h a = 0.1 (curve 1) to 0.9 (curve 2); curve 3 is determined by the regression formula for an isolated fiber, Equation 7 [32]; and curve 4 is determined by the Lamb formula for an isolated fiber [33].h a = 0.1 (1), 0.01 (2), 0.001 (3); 4 is given by the regression formula for an isolated cylinder, Equation (7) [32]; 5 is given by the Lamb formula for an isolated cylinder [33]; horizontal asymptotic lines ()-values for the Stokes flow limit at Re << 1.

Figure 3.
The dimensionless fiber drag force (F) in a row of fibers vs. the Reynolds number for different distances between the fiber axes with a/h computed with the step 0.05 starting from a/h = 0.1 (curve 1) to 0.9 (curve 2); curve 3 is determined by the regression formula for an isolated fiber, Equation 7 [32]; and curve 4 is determined by the Lamb formula for an isolated fiber [33].computed with the step 0.05 starting from h a = 0.1 (curve 1) to 0.9 (curve 2); curve 3 is determined by the regression formula for an isolated fiber, Equation 7 [32]; and curve 4 is determined by the Lamb formula for an isolated fiber [33].

Figure 3. The fiber drag force in a row of fibers vs the Reynolds number; in the case of low Reynolds numbers and sparse rows:
h a = 0.1 (1), 0.01 (2), 0.001 (3); 4 is given by the regression formula for an isolated cylinder, Equation (7) [32]; 5 is given by the Lamb formula for an isolated cylinder [33]; horizontal asymptotic lines ()-values for the Stokes flow limit at Re << 1.

Figure 4.
The fiber drag force in a row of fibers vs. the Reynolds number; in the case of low Reynolds numbers and sparse rows: a/h = 0.1 (1), 0.01 (2), 0.001 (3); 4 is given by the regression formula for an isolated cylinder, Equation (7) [32]; 5 is given by the Lamb formula for an isolated cylinder [33]; horizontal asymptotic lines (•••••)-values for the Stokes flow limit at Re 1.
Fibers 2018, 6, 90 5 of 11 Approximations for the fiber drag force in the whole considered Re range are given in the Appendix ??.Next, as expected, the decrease in a/h results in the decrease in the lateral interaction of the neighboring fibers, and they begin to behave like an isolated cylinder.For comparison, we have also plotted (in Figures 3 and 4) the curves for an isolated cylinder by the piecewise continuous approximation formula from [32]: It should be noted that the F(Re) curve for an isolated cylinder does not reach the plateau at Re 1.At the same time, the fiber drag force in the fiber row increases more slowly with Re, as compared with that of an isolated cylinder.This is explained by the mutual hydrodynamic influence of the fibers.In the next Section, the calculated flow field is used for modeling the diffusion transport in the row of fibers.

Convection-Diffusion Mass Transfer in a System of Absorbing Fibers
To determine the distribution of the solute concentration S in the flow, the elliptic equation of convection-diffusion is solved: where Pe = 2aU/D = Re(U)Sc is the fiber diameter-based diffusion Peclet Sc = ν/D is the Schmidt number, D is the solute diffusion coefficient, and u • ∇C ≡ u∂C/∂x + v∂C/∂y.Here, all the values are reduced to dimensionless form by normalizing to a, U, and the inlet concentration C 0 .The condition of full retention (absorption), C = 0, is set at the fiber surface 4. At the entrance of the cell 1, the condition of uniform concentration C = 1 is set, and at the outlet boundary 3, the zero-gradient condition, ∂C/∂x = 0.At the boundaries 2, the conditions of symmetry are used.The method of solving Equation ( 8) is based on the finite-difference scheme suggested in [34] and later used in [35].
The density of the integral diffusion flux of particles per fiber length is calculated as where r and θ are dimensionless polar coordinates and Sh is the fiber mean Sherwood number (per unit length) [36].An example of the simulated profiles of concentration is shown in Figure 5 together with the streamlines.Here, the flow asymmetry and the diffusion wake evolution with Re are highly evident.The calculated curves for η versus Re are plotted in Figure 6, where for comparison, the results of calculations by the formula of Natanson [37] (written in our notations) are provided with and without allowance for flow inertia, i.e., with F = F(Re) and F = F 0 .The drag force as a function of Re is found from the solution of the Navier-Stokes equations or is estimated by the regression formulas given by Equations ( A1) and (A4) (see Appendix, Tables A1 and   A2), while the value of 0 F for Re << 1 is given in [16].The case of dense row is illustrated in Figure 7.The formula given by Equation ( 10) was originally derived in the limit of high Peclet numbers,   Pe , and for small values of  .However, the simulations performed for a wide range of parameters reveal the wider range of validity of this formula.For estimations, it can be used at high Schmidt numbers, at high and intermediate Peclet numbers, and at h a < 0.5, i.e., for sparse and moderately dense rows of fibers.As follows from the figures, the formula of Natanson has better consistency with simulations for b = 0.2 and Sc = 1000.For this particular case, the following approximation can be used: For the range 100  Re  300 and for non-dense rows, we suggest a simple formula, which is obtained using the corresponding regression formula for the drag force in Equation ( 5): It is seen from the above figures that the slope of the   Pe  curves varies with Re (the curves become nonlinear), which is explained by mixing in the wake region (Figure 5), which leads to the increase of the diffusion flux to the fiber.The slope varies also as Pe decreases at any Sc, which is associated with a transition from a convective-diffusion transport to a fully diffusion transport when the curves reach the plateau at Pe  0 [38].In this case, under the condition 0  C at the fiber surface, the fibers completely absorb the arriving solute and the fiber retention efficiency is equal to the inlet flux, which in dimensionless form, is equal to the distance between the axes of the neighbor fibers.This case is illustrated in Figure 6a, where the curve 4 was plotted by the asymptotic formula for the fiber retention efficiency, derived at Pe << 1 in [38] for a uniform creeping flow with    The drag force as a function of Re is found from the solution of the Navier-Stokes equations or is estimated by the regression formulas given by Equations (A1) and (A4) (see Appendix ??, Tables A1  and A2), while the value of F 0 for Re 1 is given in [16].The case of a dense row is illustrated in Figure 7.The formula given by Equation ( 10) was originally derived in the limit of high Peclet numbers, Pe → ∞ , and for small values of η.However, the simulations performed for a wide range of parameters reveal the wider range of validity of this formula.For estimations, it can be used at high Schmidt numbers, at high and intermediate Peclet numbers, and at a/h < 0.5, i.e., for sparse and moderately dense rows of fibers.As follows from the figures, the formula of Natanson has better consistency with simulations for b = 0.2 and Sc = 1000.For this particular case, the following approximation can be used: η = 0.67/Pe + 2.8(F/4π) 1/3 Pe −2/3 .For the range 100 ≤ Re ≤ 300 and for non-dense rows, we suggest a simple formula, which is obtained using the corresponding regression formula for the drag force in Equation ( 5   It is seen from the above figures that the slope of the η(Pe) curves varies with Re (the curves become nonlinear), which is explained by mixing in the wake region (Figure 5), which leads to the increase of the diffusion flux to the fiber.The slope varies also as Pe decreases at any Sc, which is associated with a transition from a convective-diffusion transport to a fully diffusion transport when the curves reach the plateau at Pe → 0 [38].In this case, under the condition C = 0 at the fiber surface, the fibers completely absorb the arriving solute and the fiber retention efficiency is equal to the inlet flux, which in dimensionless form, is equal to the distance between the axes of the neighbor fibers.This case is illustrated in Figure 6a, where the curve 4 was plotted by the asymptotic formula for the fiber retention efficiency, derived at Pe 1 in [38] for a uniform creeping flow with u = (1, 0) in a sparse row of fibers: where K 0 (z) is the modified Bessel function of the imaginary argument.It should be noted that small Peclet numbers are realized in gases at low Reynolds and low Schmidt numbers.

Discussion
In this work, the hydrodynamic resistance and convection-diffusion mass transfer in the isolated rows of parallel circular fibers arranged normally to the direction of the external laminar flow, which is used as the simplest model for a hollow-fiber membrane contactor with regions of constricted flow, were simulated.The computations were preformed for the Stokes flow and inertial flow regimes for the Reynolds numbers up to Re = 300.The dimensionless fiber drag force per unit length F was calculated, which determines the dimensionless pressure drop in the fiber row system, and a correlation for the F(Re) dependence was proposed.Also, the dimensionless solute retention efficiency of a fully absorbing fiber in a row of fibers was calculated, using η = 2πSh, where Sh is the fiber Sherwood number, versus the Reynolds number for the prescribed Schmidt number values of Sc = 10, 100, and 1000.It was shown that the Natanson formula, which simultaneously relates the fiber retention efficiency η with F, Re, and Sc, and which was originally derived in the limit of high Peclet numbers for the case of small retention efficiency, is applicable for small and intermediate Reynolds numbers; intermediate and large Peclet numbers, where Pe = Re × Sc; and for sparse and moderately dense rows of fibers for a/h < 0.5.The results obtained allow one to estimate the solute retention efficiency of a streamlined hollow fibers in a contactor with simultaneous consideration of contactor parameters, such as fiber packing density, and process conditions characterized by dimensionless criteria-the Reynolds, Peclet, and Schmidt numbers.The results can find application in the field of separation and purification with absorbing fibrous media.The maximum and mean average percentage errors of the derived correlations do not exceed max (ε) = 2 and <ε> = 0.1 %, respectively.

Figure 4 .
Figure 4.The fiber drag force in a row of fibers vs the Reynolds number: in the case of intermediate Reynolds numbers, the h a values are marked on the curves.

Figure 2 .
Figure 2. The fiber drag force in a row of fibers vs. the Reynolds number: in the case of intermediate Reynolds numbers, the a/h values are marked on the curves.

Figure 2 .
Figure 2. The dimensionless fiber drag force (F) in a row of fibers vs the Reynolds number for different distances between the fiber axes with h acomputed with the step 0.05 starting from h a = 0.1 (curve 1) to 0.9 (curve 2); curve 3 is determined by the regression formula for an isolated fiber, Equation 7[32]; and curve 4 is determined by the Lamb formula for an isolated fiber[33].

Figure 3 .
Figure 3.The fiber drag force in a row of fibers vs the Reynolds number; in the case of low Reynolds numbers and sparse rows:h a = 0.1 (1), 0.01 (2), 0.001 (3); 4 is given by the regression formula for an isolated cylinder, Equation (7)[32]; 5 is given by the Lamb formula for an isolated cylinder[33]; horizontal asymptotic lines ()-values for the Stokes flow limit at Re << 1.

Fibers 2018, 6 ,Figure 2 .
Figure 2. The dimensionless fiber drag force (F) in a row of fibers vs the Reynolds number for different distances between the fiber axes with h acomputed with the step 0.05 starting from h a = 0.1 (curve 1) to 0.9 (curve 2); curve 3 is determined by the regression formula for an isolated fiber, Equation 7[32]; and curve 4 is determined by the Lamb formula for an isolated fiber[33].

Fibers 2018, 6 Figure 5 .
Figure 5.The flow and concentration C fields past a row of fibers at different values of the Reynolds number, where Sc = 1000, h a = 0.5, and the streamline inlet ordinates are 10 4 , 0.2, 0.5, 1, and 1.5.
 z K 0 is the modified Bessel function of the imaginary argument.It should be noted that small Peclet numbers are realized in gases at low Reynolds and low Schmidt numbers.

Figure 5 .
Figure 5.The flow and concentration C fields past a row of fibers at different values of the Reynolds number, where Sc = 1000, a/h = 0.5, and the streamline inlet ordinates are 10 −4 , 0.2, 0.5, 1, and 1.5.

FibersFigure 6 .F
Figure 6.Fiber retention efficiencies per unit length in a row of fibers with different interfiber distances: (a) h a = 0.2 and (b) 0.5, vs Reynolds number, found from the combined numerical solution of the Navier-Stokes and convection-diffusion equations (curves 1): Schmidt numbers are marked on the curves.Curves 2 were plotted by Equation (10) with allowance for   Re F F  ; curves 3 were plotted by Equation (10) with F = 0 F for the Stokes flow approximation at Re << 1; curves 4 were plotted by the asymptotic formula (Equation (12)) for Pe << 1.

Figure 6 .
Figure 6.Fiber retention efficiencies per unit length in a row of fibers with different interfiber distances: (a) a/h = 0.2 and (b) 0.5, vs. Reynolds number, found from the combined numerical solution of the Navier-Stokes and convection-diffusion equations (curves 1): Schmidt numbers are marked on the curves.Curves 2 were plotted by Equation (10) with allowance for F = F(Re); curves 3 were plotted by Equation (10) with F = F 0 for the Stokes flow approximation at Re 1; curves 4 were plotted by the asymptotic formula (Equation (12)) for Pe 1.

Figure 6 .
Figure 6.Fiber retention efficiencies per unit length in a row of fibers with different interfiber distances: (a) h a = 0.2 and (b) 0.5, vs Reynolds number, found from the combined numerical solution of the Navier-Stokes convection-diffusion equations (curves 1): Schmidt numbers are marked on the curves.Curves 2 were plotted by Equation (10) with allowance for   Re F F  ; curves 3 were plotted by Equation (10) with F = 0 F for the Stokes flow approximation at Re << 1; curves 4 were plotted by the asymptotic formula (Equation (12)) for Pe << 1.

Figure 7 .
Figure 7. Fiber retention efficiencies per unit length in a dense row of fibers with h a = 0.8 vs Reynolds number, found from the combined numerical solution of the Navier-Stokes and convection-diffusion equations (curves 1): Schmidt numbers are marked on the curves; curves 2 were plotted by Equation (10) with F = 0 F for the Stokes flow at Re << 1; curves 3 were plotted by the asymptotic formula (Equation (12)) for Pe << 1.

Figure 7 .
Figure 7. Fiber retention efficiencies per unit length in a dense row of fibers with a/h = 0.8 vs. Reynolds number, found from the combined numerical solution of the Navier-Stokes and convection-diffusion equations (curves 1): Schmidt numbers are marked on the curves; curves 2 were plotted by Equation (10) with F = F 0 for the Stokes flow at Re 1; curves 3 were plotted by the asymptotic formula (Equation (12)) for Pe 1.