The effect of particle size on flow in a continuous oscillatory baffled reactor using CFD

Experimental and numerical studies of flow characteristics in continuous oscillatory baffled reactors (COBR) have mainly been based on single liquid phase in the past decades. This work, for the first time, investigates the effects of particle size on axial dispersion and evaluates the residence time and velocity experienced by solid particles in a COBR, along with their impact on solid suspension. Computational fluid dynamics (CFD) methodology was employed for this work by coupling a primary Eulerian liquid phase with a secondary discrete Lagrangian phase consisting of solid particles of given density and size. The data not only provide insights on how particles behave in a COBR, but also assist the design and development of COBR.


Introduction
Oscillatory Baffled Reactors (OBR) offer uniform mixing 1 and linear scale up 2 , making them an attractive alternative to stirred tank reactors for research and industrial applications in reaction [3][4][5][6][7] and crystallization processes [8][9][10][11][12][13][14][15][16] . The characterization of flows in COBR has extensively been reported using single phase [17][18][19][20][21][22][23][24][25][26][27][28][29] in the past decades. Recent work by Ejim et al. 30 highlighted the differences and the knowledge gap on the design of COBRs for multiphase flow processes using correlations obtained from single phase studies; this was further emphasised by Kacker et al. 31 who reported that not only the optimal operating conditions for minimal axial dispersion involving solids were different from that of single phase, but also longer times were spent by solids in a COBR, highlighting the need to properly address the effect of different solid particles on axial dispersion. Baptista et al. 32 analyzed the behavior of suspended solid particles of different sizes and densities in a baffled reactor; however, their system did not include oscillatory flow and their findings were inconclusive, as the interaction among particles was too significant for the effects of size and density to be evaluated.
From a modelling standpoint, Computational Fluid Dynamics (CFD) solvers have been developed and used to simulate hydrodynamic flow profiles in OBRs/COBRs for multiple purposes [33][34][35][36][37][38][39] . Again, all these studies were performed under a single liquid phase framework, including the work by Mazubert et al. 40 , who made use of discrete particle tracking of a secondary phase to measure concentration profiles and analyze the performance of different geometric designs; however, their secondary phase consisted of massless particles that essentially followed the velocity field of the continuous Eulerian phase.
The present CFD work involves two phases and aims to investigate the effect of particle size on axial dispersion in COBRs, while evaluating the residence times and velocities experienced by these solid particles and their impact on solids suspension. To the authors' best knowledge, this is the first study of its kind in the area of COBR where a continuous Eulerian phase is coupled with a discrete Lagrangian (solid) phase.

Problem definition
The dimensionless numbers that govern the conditions of the flow in a COBR are the net flow Reynolds numbers (Ren = unet-inletρD/μ), the oscillatory Reynolds number (Reo = ωxoρD/μ), the Strouhal number (St = D/4πxo), the ratio of the area of the orifice over the area of the tube, known as the restriction ratio (α = Db 2 / D 2 ), and the velocity ratio (Ѱ = Reo/Ren); where ρ and μ are, respectively, the fluid density (kg m -3 ) and dynamic viscosity (kg m -1 s -1 ), D is the diameter of the tube (m), Db the diameter of the baffle hole (m), ω=2πf the oscillation angular frequency (rad s -1 ), xo the oscillation center-to-peak amplitude (m), f the oscillation frequency (Hz) and unet-inlet is the net inlet velocity (m s -1 ).
It should be noted that unet-inlet was carefully labeled to avoid any confusion with unet as unet ≠ unet-inlet. The net volumetric flow rate (Q) is constant and the net velocity at the inlet (unet-inlet) is thus calculated as Q/A. However, while both Q and unet-inlet are constants, the net velocity changes along the length of the reactor due to its smooth-edged baffles, as shown in Figure 1, e.g. the volume of fluid flowing through a baffled constriction is 25% less than that of a cylinder of a diameter D, indicating V ≠ L•A in this case; subsequently, the velocity through orifices (unet-baffle) is defined as Q/Ab and the mean net velocity (unet) of the system is within the range Q/A ≤ unet ≤ Q/Ab. If the volume of the device is known, unet is calculated as QL/V.
The target device of this study is a NiTech DN15 COBR reactor, as shown in Figure 1, the geometric dimensions of the DN15 and all design details were provided by the manufacturer, Alconbury Weston Ltd (http://www.a-w-l.co.uk/); the total length of the reactor is 752 mm, containing 32 baffle-cells. The operating oscillatory conditions were chosen partially based on the characteristics of the simulated particles, the simulated domain and literature work.
Kacker et al. 31 used a net flow rate of 100 ml min -1 and identified an oscillatory amplitude of 2 mm and frequency of 2 Hz as optimal conditions for solid suspension and near plug flow behaviour of melamine crystals (mean particle size = 100 μm). Hence, Q = 100 ml min -1 and f = 2 Hz were selected for this study.
Extensive literature has reported a proportional relationship between oscillatory amplitude and axial dispersion [19][20][21]27,31,41 . However, while minimal dispersion and near plug flow behaviour are desirable, Oliva et al. 42 stated that the minimal energy required to ensure solid suspension should be considered. For this reason and considering paracetamol particles of up to 150 μm in diameter were simulated in this work, a moderate oscillatory amplitude of xo = Db (7 mm) was selected based on the work by González-Juárez et al. 39 . It should be noted that while the chosen 7 mm amplitude is within the reported range for RTD studies in DN15 31,42 , it is lower than that used in crystallization processes in COBRs where it commonly ranged from 12 to 30 mm 15,16,43-45 [references].
Due to the nature of oscillatory flow in COBRs, forward and backward mixing are generated during oscillation, resulting in particles flowing in and out of a control domain. It is thus crucial to select the injection point for solid phase as well as the measuring points to ensure that open boundary effects on particles, due to the inlet and outlet, are minimized. In terms of the injection point, it was set at the baffle-cell number 15, i.e. 352.5 mm from the inlet, which ensured that less than 0.1% of the injected particles left the system through the inlet for all the simulated conditions, see Table 1.

Determination of axial dispersion
In a tubular reactor, mixing is commonly quantified by axial dispersion coefficient (Da), which describes the degree of spreading (in the axial direction) of a tracer injected upstream as a pulse (ideally). Analogous to the molecular diffusion model, Levenspiel and Smith 46 proposed the following equation to evaluate the axial dispersion coefficient: where C is the tracer concentration as a function of time (t) and position (x) and U is the mean net flow velocity of the system (U = QL/V). Equation (1) was originally derived for a single phase flow and could also be used for two phase (solid -liquid) cases. When the concentration of a liquid tracer is defined as CL = mL/VL and the concentration of solids as CS = mS/(VL + VS), eq. (1) becomes independent of the volume of the secondary solid phase (VL + VS ≈ VL) when VS << VL, which is the case in the present study. This is also consistent with literature, e.g. the work of Ejim et al. 30 and Kacker et al. 31 . If a perfect input pulse injection is assumed, typical boundary conditions for eq. (1) are: where n is the volume of tracer/secondary phase injected, A the cross-sectional area of the device and δ(x) a Dirac delta function. Thus, the analytical solution to equations (1) and (2) at fixed values of Da and U is given by Under the assumption of a perfect pulse injection, the plug flow with axial dispersion model can be re-derived and solved based on an inverse Peclet number (Pe = UL/Da) as: where σθ 2 is the dimensionless variance, σ 2 the variance and ( ) t the mean residence time of the tracer concentration as: where N is the total number of measured concentration data. Eq. (4) is used to evaluate axial dispersion coefficient under the perfect pulse method (PPM).
A perfect input pulse is virtually unachievable. Hence, an imperfect pulse method (IPM) was proposed 47 where the concentration profile of the tracer/secondary phase is measured at two downstream points of the tracer injection, i.e. C1(t) and C2(t); thus the form of the impulse becomes irrelevant. This method was firstly implemented in OBRs by Mackley and Ni 1 , who adopted the solution of Göeble et al. 48 and used a normalized concentration E(t) for better comparison among results: Mackley and Ni 1 suggested that the normalized concentration measured at an upstream point where the distance from the injection point (L) is essentially the distance between measuring points (1) and (2). The normalized concentration predicted at point (2), E2′(t), is compared with the measured normalized concentration at such point, E2(t), and the axial dispersion coefficient is fitted in order to satisfy the target function: where N is the total number of normalized concentration data. The optimal axial dispersion coefficient is obtained when the target function (10) is minimized. While the value of the mean net flow velocity (U), as aforementioned, can be assumed as U = unet, a better method is to calculate the time it takes for the tracer/secondary phase to travel from measuring point (1) to (2) as: The imperfect pulse method 20

Computational simulation set-up
All numerical simulations were performed using ANSYS® Fluent 16.0 CFD package, which discretizes a computational domain using finite volume methodology in order to solve the flow field of a continuous phase. Additionally, Fluent allows for Lagrangian particle tracking by implementing a so-called Discrete Phase Model (DPM) as an add-on to an existing Eulerian phase, this capability was utilized to model massless (tracer) and solid particles.

Numerical model for Eulerian phase
The fluid selected for this study was water (ρ = 998.2 kg m -3 , μ = 1.003•10 -3 kg m -1 s -1 ); timedependent incompressible three-dimensional Navier-Stokes equations were solved as: All simulations were performed utilizing a pressure-based segregated solver along with the SIMPLE pressure-velocity coupling algorithm. Spatial discretization of the momentum equation was performed using a second order upwind scheme; pressure at faces of the grid was interpolated using a second order scheme and time was discretized using a first order implicit scheme. The time-step was set to 2 ms, ensuring a good number of time-steps per oscillatory cycle of 250, which is higher than the norm reported in literature 39,40 . The average and maximum values of the Courant-Friedrichs-Lewy (CFL) coefficient were consistently kept below 2 and 20 respectively.
Because of the low net flow and oscillatory Reynolds numbers of the conditions simulated, a laminar solver was selected. This is also in agreement with existing literature 35,[38][39][40][52][53][54][55] . The impact of inlet boundary conditions on the main flow was minimized by imposing a fully developed parabolic profile: . A constant gauge pressure of 0 Pa was set for the outlet boundary; operating conditions were set at 300 K and 101325 Pa.

Mesh sensitivity test
A 5-baffle-cell tube geometry was employed for the performance of a mesh sensitivity study, as illustrated in Figure 3. This analysis was previously undertaken during a power density study in COBRs 56 [UPDATE], under its most adverse conditions (Q = 50ml min -1 , f = 8Hz, xo = 14mm). The time-step was set to 0.5 ms for an oscillatory cycle consisting of 250 timesteps. Considering global mesh refinement, five meshes of different resolutions were compared and each simulation was run for 24 oscillatory cycles. The variables compared between meshes were pressure drop vs time profiles (∆p(t) = p1(t) -p2(t)) and velocity magnitude vs time profiles extracted at lines 1 & 2 and planes 1 & 2 as shown in Figure 3.
These profiles were cycle-averaged and the resultant profiles (of the duration of an oscillatory cycle) were compared using the coefficient of determination (R 2 ): where SSres is the sum of squares of residuals between the target profile (that from mesh #1) and the profile under evaluation (from mesh #j) and SStot is the total sum of squares of the target profile. Subscripts i and n represent a single data point and the total number of data points of a certain profile respectively, while j is the index of a certain mesh and ɸ the property under evaluation. The profiles extracted from each mesh were compared to those from mesh #1. The results of this mesh sensitivity analysis are summarized in Table 2; Mesh #2 (in bold) was selected on the balance of accuracy and efficiency, and its density is above the norm reported in literature 36,38,52,54,[57][58][59] . All meshes were created on ANSYS ICEM containing only hexahedral elements and were O-grid structured.

Numerical model for Lagrangian phase
The discrete solid phase was mono-sized spherical paracetamol (ρ = 1263 kg m -3 ) particles of diameter (Dp) 50, 100 and 150 μm, while liquid phase information was obtained from discrete massless particles that act as a perfect tracer as they move according to the flow field of the continuous liquid phase. The trajectory of each discrete particle is predicted by integrating the force balance on the particle as: where mp, p u  and p ρ are, respectively, the mass, velocity and density of the particle. The second term in the right-hand side of eq. (15) accounts for the force due to the weight of the particle and the buoyancy effect, the first term D F  is the drag force defined as: where Ap is the cross-sectional area of the particle and CD is the drag force coefficient, calculated as the spherical drag law proposed by Morsi  , the former accounts for the force required to accelerate the fluid surrounding the particle and the latter is the resultant force from the pressure gradient along the fluid flow around the particle: where D Dt is the material derivative. The position of each particle ( ) Equations (15) and (19) are integrated using a trapezoidal discretization scheme with the same time-step as the Eulerian phase (2 ms). All particles were modelled as perfect spheres and were released at a cross-sectional plane at the middle of a selected baffle-cell; this is the so-called "surface injection". In order to cope with the potential computational limitation of modelling too many particles, Fluent tracks so-called "parcels". A parcel may contain multiple particles; its position is defined by a tracked representative particle and its diameter is that of a sphere whose volume is the ratio of the total parcel mass to particle density.
However, in this work, in order to model and predict the behavior of individual particles, the mass of each parcel was set as that of a single particle, i.e. each parcel contained one particle and thus the concept of parcel and particle are interchangeable in this study. The number of particles released in the system was set to 4050; further analysis on the sensitivity of the number of tracked particles on results will be discussed on section 5.1. No particle -particle interaction or particle diffusion in the liquid phase were included in the model.

Effect on number of simulated particles
Discrete particles were injected into a cross-sectional plane in the middle of a baffle-cell containing 4050 computational cells; the "surface" injection type releases one particle per computational cell, thus 4050 particles were injected. This number is significantly larger than what is reported in literature, e.g. Mazubert et al. 40 examined the effect of the number of injected particles (2484 and 4968) on the simulated results and found that the difference between the two was negligible. Taking a conservative approach, particle numbers of 4050 and 8100 were examined and compared in this work; for the latter, the "surface" injection was simultaneously performed at two cross-sectional planes in the middle of the 15 th bafflecell, with a distance of 0.75 mm apart. This analysis was performed for massless particles under the operating conditions of run #2 (see Table 1). Figure   RTD curves measured at three different baffle-cells during run #2 are presented in Figure 5 (left).

Figure 5. RTD curves measured at different baffle-cells at operating conditions of run #2 (left) and RTD area under the curve (Co) with length for all simulated conditions (right)
A constant Co value for each simulated run is a good guide for selecting measuring points, see the framed square in Figure 5 (right). Consequently, RTD data obtained from baffle-cells (17) to (27) (0.047 -0.282 m from the injection point) was selected for analysis, while the remaining concentration profiles measured at baffle-cells (28) to (32) were discarded as the effect of the open boundary was too large for reliable C(t) curves to be measured.  Figure 6 shows the axial dispersion coefficients obtained from both the PPM and IPM methods as a function of length for all runs (refer to  (20) to (27) were set as C2.

Perfect vs Imperfect pulse method
Although Da values computed via IPM fluctuate around the asymptotic value, they are much more stable than their counterparts and are thus chosen as the final results in this study.

Figure 6. Da calculated from RTD curves measured at different lengths (from the injection point) using the imperfect (IPM) and the perfect pulse (PPM) methods for all runs simulated.
During this analysis, it was observed how the residual errors from target function (10) were consistently higher when U was assumed to be equal to QL/V, as opposed to when it was calculated with eq. (11). Figure 7 shows the velocity values obtained during the fitting process of IPM, calculated via eq. (11) at different lengths of the reactor. In reality, U from eq. (11), presented in Figure 7, represents the mean net velocity particles experienced while travelling from measuring point (1) to point (2). This velocity becomes more smoothed out as measuring point (2) is moved along the length of the device at a fixed measuring point (1), thereby increasing the length over which RTD curves are examined. Our preliminary results clearly indicate that while U = QL/V is a fair estimation of liquid phase velocity, the velocity of a secondary solid phase is dependent on particle size, with velocity of small particles (50 μm diameter) close to that of the liquid phase and that of larger particles being significantly smaller. As a result, the assumption of U = QL/V when applying the IPM was disregarded.
Note that the discussion above refers to a secondary phase containing mono-size particles similar to those simulated in this work.

Figure 7.
Velocity calculated with eq. (11) with measuring point (1) fix at baffle-cell (17) and measuring point ranging from baffle-cell (19) to (27)      The results they obtained from injection port 1 and measuring port 1 (L ≈ 0.72 m) were chosen as the basis for comparison; the profile measured at baffle-cell (27) of our device (L = 0.282 m) shown in Figure 8 (right) has a remarkable similarity both in magnitude and shape, hence validating the simulated results of this work.

Effect of size of particle on axial dispersion and residence time
The effect of particle size on axial dispersion is graphically presented in Figure 9 (left), displaying profiles of E(θ) vs θ of the baffle-cell (27) for all runs performed (see Table 1).
While the impact of particle size on axial dispersion of the secondary phase is qualitatively minimal (Figure 9 left), the asymptotic Da values from Figure (6) decrease with the increase of particle size ( Table 3). The mean residence time required for particles to reach baffle-cell  (Table 3).

Effect of particle size on their velocity and suspension
In this work, the time-dependent velocity magnitude (um) and the velocity in the axial direction (ux) were calculated as: where the index i represents a specific particle. Because the total number of particles is a µm respectively are clearly seen as particle size increases.   The minimum transport velocity required for the suspension of slurry in a horizontal tube (umin-h) is given by the modified Durand equation 61,62 : where D is the diameter of the tube, dp the diameter of the particle and Cmh an empirical constant that ranges from 0.4 to 1.5; taking a conservative approach, a Cmh of 1.5 was used in our calculations. Table 4 shows the umin-h required by each set of simulated solid particles, along with the percentage of the time that the oscillatory inlet velocity is higher than the minimum transport velocity, i.e.  The minimum transport velocity required for slurry suspension in a horizontal tube increases by 20% as particles grow from 50 µm to 150 µm of diameter. The inlet oscillatory velocity is greater in magnitude than umin-h more than 50% of time for all simulated particles; hence, it is expected for all particles to stay suspended throughout their journey downstream the reactor.
However, results point to a decrement in the degree of suspension as particles grow in size. As expected, Figure 14 shows that the average position in the y-axis decreases as particle size increases. While the suspension for liquid phase and small solid particles (50 µm) occurs at the centre of the device (y = 0), particles of 100 and 150 µm lose height and stay somewhat closer to the bottom wall of the reactor (y = -7.5 mm) as they move axially downstream; these bigger particles are a good example of heterogeneous suspension 61 . Complete settlement of particles was not seen even for the largest particles injected.

Conclusions
By coupling a primary Eulerian liquid phase with a secondary discrete Lagrangian phase, we have, for the first time, reported a detailed analysis on the effect of particle size on axial dispersion, residence times and velocities experienced by solid particles in a COBR, as well as their impact on solid suspension. Results show a decreasing trend in oscillatory axial velocity as particle size increases, leading to smaller axial dispersion and longer residence times. These findings agree with the work by Ejim et al. and Kacker et al. 30,31 .
On the determination of axial dispersion of secondary phase, two methodologies were utilized in this study: the perfect pulse method and the imperfect pulse method. The latter provided constant results at different lengths of the device for all the simulated cases, while the former did not. This is most certainly due to the formulation of the IPM that avoids the assumption of a perfect pulse injection of the secondary phase.