The effect of Stefan flow on the drag coefficient of spherical particles in a gas flow

Particle laden flows with reactive particles are common in industrial applications. Chemical reactions inside the particle can generate a Stefan flow that affects heat, mass and momentum transfer between the particle and the bulk flow. This study aims at investigating the effect of Stefan flow on the drag coefficient of a spherical particle immersed in a uniform flow under isothermal conditions. Fully resolved simulations were carried out for particle Reynolds numbers ranging from 0.2 to 14 and Stefan flow Reynolds numbers from ( −1 ) to 3, using the immersed boundary method for treating fluid-solid interactions. Results showed that the drag coefficient decreased with an increase of the outward Stefan flow. The main reason was the change in viscous force by the expansion of the boundary layer surrounding the particle. A simple model was developed based on this physical interpretation. With only one fitting parameter, the performance of the model to describe the simulation data were comparable to previous empirical models. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ ) c n c l t a n t s s c l t s i b i


Introduction
Many industrial applications involve particle laden flows with reactive particles, such as combustion of solid fuels, catalytic cracking and drying applications. Unlike ordinary particle-laden flows, reacting particles exchange mass with the surrounding fluid. A Stefan flow, induced by chemical reactions inside or at the surface of the particle, has effects on the gas-solid interaction, i.e. momentum ( C D -drag coefficient), heat ( Nu -Nusselt number) and mass transfer ( Sh -Sherwood number) between the particle and the bulk flow ( Hayhurst, 20 0 0;Yu and Zhang, 2009;Yu et al., 2013;Kalinchak, 2001 ). This can be exemplified by gasification and combustion processes, where, upon being released into the hot environment, fuel particles undergo fast devolatilization that results in a pronounced gas stream leaving the particles. Although momentum, heat, and mass transfer could be affected by the Stefan flow, as a first step, we focus on the effect of Stefan flow on C D in isolation from the effects of heat and mass transfer in this study. * Corresponding author.
Resolved simulations of multiphase reactive flows demand high computational resources due to its complexity and the multi-scale nature of the processes. The smallest scale in such systems typically corresponds to the scale of the particles and their boundary layers ( 10 −6 − 10 −3 m), while the largest scales are set by the entire reactor, which typically contains millions of reactive particles and has a length scale ( 10 0 − 10 2 m) that is several orders of magnitudes larger than the particle scale. Therefore, it is impractical to carry out particle resolved simulations for a large domain. Instead, it is useful to develop constitutive models based on the results from particle resolved simulations of single or multiple particles, which can then be implemented in large scale reactor simulations that do not resolve the individual particles. Compared with the many particle-resolved simulations in the literature, only a few studies have used their results to develop models suitable to use in large scale simulations (e.g. models for Stefan flow developed by Miller and Bellan (1999) and Kestel (2016) , while models taking into account particle porosity and particle shape are presented in Wittig et al. (2017) and Richter and Nikrityuk (2012) , respectively.) Previous studies on Stefan flow effects mainly investigated droplet evaporation/condensation ( Bagchi et al., 2001;Renksizbulut and Yuen, 1983;Dukowicz, 1984 ) and suction/blowing effects ( Chuchottaworn et al., 1983;Dukowicz, 1982;Cliffe and Lever, 1985 ). Models developed for the drag coefficient of evaporating/condensing droplets are based on both experimental and simulation data. Recently, the performance of the model by Renksizbulut and Yuen (1983) was assessed for a char particle during oxy-fuel combustion ( Farazi et al., 2016 ). The model contains a case-specific blowing number and had to be adjusted by introducing a new blowing number. However, some studies have proposed more general models for the drag coefficients of a reacting particle, based on the suction/blowing effect directly. In early models, the mass flux inward/outward (hereafter called Stefan flow) was represented by 'a non-dimensional blowing number ( )', which is the ratio of Stefan flow velocity and slip velocity ( = U s f /U ∞ ) ( Cliffe and Lever, 1985 ). More recently, the Stefan velocity has been non-dimensionalized by the Stefan Reynolds number, Re sf , which is based on particle radius ( R ), Stefan velocity ( U sf ) and fluid viscosity ( ν) ( Kestel, 2016 ): (1) Another relevant Reynolds number is the particle Reynolds number, Re , which is based on the particle slip velocity ( U ∞ ), such that U s f,r = Re s f /Re . Dukowicz (1982) developed an analytical relation for the drag of a spherical solid particle with suction/blowing in creeping flows ( Re → 0). For higher Re , a number of works addressed the effects of Stefan flow on the drag coefficient ( Cliffe and Lever, 1985;Miller and Bellan, 1999;Kestel, 2016;Nour et al., 2017 ). Miller and Bellan (1999) developed an empirical model based on the numerical simulation results of Cliffe and Lever (1985) for an isothermal flow around a sphere. Kurose et al. (2003) has modified the model coefficients of the same model to fit the data for an outflow in a linear shear flow around a solid sphere. Later, another empirical model was introduced by Kestel (2016) , which is applicable for the wider range of mass fluxes that appeared in a 200 MW commercial gasifier data. It is apparent that the change of drag coefficient due to Stefan flow cannot be neglected. However, available models are not based on physical observations, and they rely on a number of fitting parameters. In addition, none of the models are suitable for negative Stefan flows (suction).
This study investigates the interaction between a gas flow and an embedded reacting particle that experience a Stefan flow. The main aim is to develop a physics-based simple model describing the change of the drag coefficient due to the Stefan flow for a particle in an isothermal flow. Direct numerical simulations that resolve the boundary layer at the particle surface were carried out for a laminar flow surrounding a stationary particle with either an outgoing or an incoming Stefan flow. Simulation results were analyzed and a model was developed with a physical interpretation from the simulations. The developed model and two previous models from the literature ( Miller and Bellan, 1999;Kestel, 2016 ) were compared with the simulation results. The range of particle Reynolds numbers ( Re ) in this study is limited to the conditions relevant to entrained-flow gasification or pulverized combustion.

Methodology
The numerical simulations considered a static particle in a uniform isothermal flow. The generation and consumption of gas in the solid phase were considered as a uniform outgoing or incoming mass flux at the particle surface in the surface-normal direction. In all of the simulations performed here, the Reynolds number is smaller than the critical Reynolds number that yields von Karman oscillations. This means that there are no transients in the flow, and hence, a steady state solver can be used.

Governing equations
Steady state simulations were carried out under isothermal conditions, with the gas phase assumed to be incompressible. The discrete phase was described as a static spherical particle with constant size. The gas phase is governed by mass conservation, and momentum conservation, where ρ is the density of the fluid, − → u is velocity vector, p is pressure and μ is dynamic viscosity. Eqs. (3) and (4) were discretized with the finite volume method using second-order schemes.

Boundary conditions
The slip velocity between the particle and the bulk gas was set as the inlet velocity at the front boundary (left side of the calculation domain in Fig. 1 ). An 'outflow' boundary condition (i.e. zero velocity gradient) was applied at the back boundary (right side of the calculation domain in Fig. 1 ). The side boundaries of the domain were treated as 'slip walls'. A 'slip wall' boundary condition enforces both the velocity component normal to the wall and the gradients of the other velocity components in the normal direction to be zero. Boundaries along the symmetry axes were considered as 'symmetric' boundaries, which means that the component of velocity normal to the symmetry plane is zero and that the gradient of all the other properties normal to the plane is zero.
The immersed boundary method (IBM) was used at the surface of the particle. The current work used the discrete forcing approach ( Mittal and Iaccarino, 2005 ), which uses the direct imposition of boundary conditions ( Jasak et al., 2014 ), and the presence of the immersed surface/body is formulated through the boundary conditions. The value of any parameter inside the cells that contain the immersed boundary was calculated by interpolating values at the immersed boundary points and the neighbour cells ( Fadlun et al.,20 0 0 ). To implement Stefan flow, the velocity is fixed (Dirichlet boundary condition) at the immersed boundary normal to the particle surface as: where integration is over the surface S of the particle, − → n is unit vector in the direction normal to the surface element dS and ˙ m is mass flow rate due to the Stefan flow. Furthermore, for pressure the gradient is set to zero at the immersed boundary (Neumann boundary condition). The treatment of Dirichlet and Neumann boundary conditions for an immersed boundary method in foam-extend is shown in the Appendix A ( Jasak et al., 2014 ).

Calculation conditions and procedure
In this work, we used the OpenFoam environment, called foam-extend -3.2 ( Weller et al., 1998 ). The numerical simulations were carried out using the incompressible, steady-state, immersed boundary solver. The solver uses quadratic interpolation ( Jasak et al., 2014 ) for the reconstruction of the solid phase boundary conditions into the closest fluid cells.
Flow conditions were selected based on practical applications of pulverized combustion and gasification at atmospheric pressure. Four different Re were selected by considering particle size (0.1-1.0 mm), slip velocity (0.5-3 m s −1 ), and gas properties of N 2 at Computational domain for the simulations, with D denoting the particle diameter, and i , i = 1 to 5 representing the coarsest mesh to finest mesh. D −x,i is the distance from the centre of the sphere to negative x-direction and D + x,i is the distance from the centre of the sphere to positive x-direction in level i (See the Table 1 ).

Table 1
Distance from the centre of the particle in diameters ( D ) in the computational domain (See Fig. 1 ).  ( Kreitzberg et al., 2016;Umeki et al., 2012 ). Since the Re was always less than 20 in this study, the flow is steady, axisymmetric and topologically similar ( Johnson and Patel, 1999 ). Therefore, only a quarter of the domain was simulated with symmetric boundaries. Initially, the domain size and mesh resolution was selected based on previous studies ( Constant et al., 2017;Richter and Nikrityuk, 2012 ) for flow around a sphere. Then, mesh refinement tests were carried out for the highest Re . Based on these tests, we arrived at five levels of refinement that were eventually used for the simulations, with the mesh size of the finest refinement being 0.01 D (see the Fig. 1 and Tab. 1 ). After the mesh refinement test, domain size tests were carried out for the smallest Re and the highest Stefan flow velocity, i.e. because the boundary layer is expected to be the largest under such condition. Based on the results (see Table 2 ), mesh 2 was selected considering accuracy and econ-  Fig. 1 , consisting of around 9.6 million cells in total. For isothermal conditions, the drag coefficient of a particle with no Stefan flow should depend only on Re . As preliminary tests, we confirmed this with two different sets of particle diameters and slip velocities at the same Re .

Estimation of the drag coefficient
The drag coefficient can be calculated as when the pressure and viscous forces are given as and respectively. Here, the integration is over the surface S of the particle. In the above, P sur and P ref are the interpolated pressure at the particle surface and in the far field, respectively, and − → n is the unit vector in the surface-normal direction. Only the components − → F P and − → F v isc in the direction of the mean flow were accounted for when calculating the drag coefficient, since the other components are canceled out due to symmetry.

Validation
The numerical implementation was validated for the estimated drag coefficient using four Re w ithout Stefan flow. The obtained drag coefficient was compared to the empirical formula of Haider and Levenspiel (1989) ,  which was derived from 408 experimental data points. Fig. 2 shows that the drag coefficients obtained from our simulations (symbols) are in agreement with this empirical formula (solid line). The data is also listed in Table B.1 . The velocity profile surrounding the particle generated by Stefan flow was validated in a quiescent fluid by comparing it to the analytical solution, where − → u d is the velocity vector at a distance d from the centre of the sphere, and − → u s f is the Stefan flow velocity vector at its surface. Fig. 3 shows the normalized drag coefficient, C D,r = C D,s f /C D, 0 , plotted against the normalized Stefan flow velocity, U s f,r = U s f /U ∞ , for different Re . Here, C D ,0 and C D,sf refer to the drag coefficients without and with Stefan flow, respectively, while U ∞ is the inlet velocity. The results show a nearly linear relationships between C D,r and U sf,r for every given Re , with the slope of the relationship getting steeper with increasing Re . According to Fig. 3 , the normalized drag coefficient was as low as 0.7 (for Re = 2.32 and U s f,r = 1 . 3 ), and is expected to decrease  even further at higher Stefan velocity. This significant reduction in drag shows the relevance of the Stefan flow in entrained flow gasification and combustion applications. Fig. 4 explores the effect of Stefan flow in more detail by showing the pressure and viscous forces separately. In all cases studied here, both with and without Stefan flow, we found that the viscous force was larger than the pressure force by a factor of roughly two, as is expected for low Re . We do see, however, that this factor is decreasing for increasing Re sf , and for much larger values of Re sf it can not be excluded that it may even be less than one. The bottom line is that a positive Stefan flow give a significant reduction of the viscous force while the pressure force remains almost constant.

Effects of Stefan flow on drag coefficient
To elucidate the observed effects, the pressure force component in the mean flow direction, (P sur f − P re f ) − → n x is shown in Fig. 5 as a function of surface angle from the front of the particle (See schematic in the inset of Fig. 1 ). The Stefan flow velocity at the surface is given as U sf / U ∞ , where the positive values indicate outgoing flows. The figure confirms the observation from Fig. 4 , i.e.; the pressure force is hardly affected by the Stefan flow and it is almost constant for a given Re . On the contrary, Fig. 4 showed that the viscous force decreased with an outgoing Stefan flow and increased with an incoming Stefan flow. To explore this effect, Fig. 6 shows the viscous stress com- ponent in flow direction as a function of the surface angle from the front of the particle. It shows that the viscous stress is actually higher at the front of the particle for the simulations with outgoing Stefan flow. On the other hand, the viscous stress behind the particle is smaller with outgoing Stefan flow. The changes in the viscous stress at the front and the back of the particle cancel each other and have no significant net effect. The shear stress at the side of the particle (40 < θ < 140) is smaller with outgoing Stefan flow. As a result, the overall viscous stress for the particle decreased under the influence of outward Stefan flow. The main factor that affects the viscous force is the velocity gradient as shown in Eq. (8) . The observation in Fig. 6 implies that the change in the boundary layer thickness is more significant than the change in velocity difference that appear at the front and back of the particle. Fig. 7 shows the flow field (i.e. velocity magnitude) with (lower half panel) and without (upper half panel) outward Stefan flow. Comparison of the flow fields showed that the boundary layer thickness increased with outgoing Stefan flow. On the contrary, the boundary layer thickness decreased with incoming Stefan flow. This change in boundary layer thickness due to the Stefan flow is clearly seen by inspecting the velocity magnitude along the y -axis crossing the centre of the sphere, as shown in Fig. 8 . For an outgoing Stefan flow (red dashed line in Fig. 8 ) we observed a slower relaxation of the velocity magnitude to the free stream velocity, while vice versa, a faster relaxation was observed for incoming Stefan flow (green dashed line). This effect can be understood as the boundary layer being pushed away from the particle surface in case of an outward Stefan flow, while it was pulled towards the surface for an inward Stefan flow. This change in boundary layer thickness with the Stefan flow affects the velocity gradient, and hence it explains the observed change in the viscous force and, consequently, also the drag.

A model for the drag coefficient with Stefan flow
In this section, a simple expression is suggested for the drag coefficient under the influence of Stefan flow for small Re . The net drag on a particle is entirely determined by the boundary layer around the particle. Any change to the boundary layer due to the presence of Stefan flow would therefore have an effect on the drag. Indeed, we observed in the previous section that Stefan flow reduced/enhanced the drag coefficient by modifying the viscous force through the expansion/contraction of the boundary layer. As a first approximation, we assume that the change in the normalized drag coefficient depends on the change in the volume of the boundary layer. By assuming that the volume of boundary layer with Stefan flow simply becomes the sum of its original volume ( V B ) and the volume of Stefan flow ( V sf ), the normalized drag coefficient can be expressed as Here, the volume of the Stefan flow is defined as the volume of fluid emitted from the particle during a typical flow time-scale, τ f , such that where the flow time-scale is given by where δ is the boundary layer thickness. We can assume δ R at small Re . Hence, Based on the above, the volumes of the Stefan flow and its approximation at low Re are now given by Furthermore, the volume of the boundary layer is given as and when δ R , Adopting the result from classical boundary layer theory, the boundary layer thickness is given by where A is a constant with a value of the order of one. Combining Eqs. (15) and (17) with Eq. (11) yields Eq. 19 is based on the assumption that Re is small enough to satisfy δ R , and it is not applicable for higher Re .  Without the assumptions of δ R , i.e. keeping the particle radius when estimating the boundary layer volumes using Eqs. (15) and (16) , the normalized drag coefficient based on Eq. (11) follows as : where A = 3 . 01 ± 0 . 13 as obtained by fitting to the simulation data using the least squares methods. The performance of the model was compared against the previous models by Miller and Bellan (1999) ( Eq. (22) ) and Kestel (2016) ( Eq. (23) ); the former reads as: which is valid for 0 ≤ Re ≤ 100 and 0 ≤ Re sf ≤ 10 ( Miller and Bellan, 1999 ) and Kestel (2016) model reads as; where a = ( 1 . 063 1+0 . 223 Re ) 0 . 568 , which is valid for 0 ≤ Re ≤ 200 and 0 ≤ Re sf ≤ 20. Fig. 10 compares the performances of three models with the data from the simulations. All the models are in good agreement with the simulation results for positive Re . The maximum error of the current model was less than 6% in the simulated range that is 0 ≤ Re ≤ 14 and (−1) Re s f 3 . However, there are two major differences between the current and previous models. First, the previous models contain several fitting parameters without clear physical background. The current model, however, contains only one fitting parameter, which is related to the relationship between Re and the boundary layer thickness ( Eq. (18) ). Moreover, the previous models by Miller and Bellan (1999) and Kestel (2016) are not applicable to negative Re sf while the current model expands to negative Re sf and is in good agreement with simulation data, at least down to Re s f = (−1) . For strongly negative Re sf , C D,r given by Eq. (20) di- verges. However for Re = 0 . 232 , Re sf has to become as small as (-7) before C D,r diverges.

Conclusions
Fully resolved numerical simulations of flow surrounding a gasemitting particle were carried out to elucidate the effect of Stefan flow on the drag acting on a particle in a uniform flow. The application of this study is limited to steady, axisymmetric flow ( Re < 14), and low Stefan flow velocity ( −1 Re s f 3 ). Results showed that the drag coefficient has a nearly linear relationship with the Stefan flow velocity. An outward Stefan flow lead to a reduction of the drag coefficient, whereas the magnitude of the reduction increases with increasing Re . For the Reynolds numbers in this study, the main reason for the reduction/increase in the drag coefficient was the change in viscous force. This was caused by the expansion/contraction of the boundary layer surrounding the particle, rather than the change in relative velocity at the particle surface.
A simple model was developed based on the change in the volume of the boundary layer due to Stefan flow. Although the model contains only one fitting parameter, it showed as good agreement with the simulation data as previous models with several fitting parameters. The proposed model also showed good agreement with the simulation data for negative Re sf while previous models could not be computed because of non-integer indexes for Re sf . Similar studies for Nusselt number and Sherwood number would be important for future works. Computing (SNIC) at High Performance Computing Center North (HPC2N). The authors thank all the staff of HPC2N for the technical assistances. Furthermore, M.U.B. thanks the Swedish for Gasification Center and its industrial and academic partners for financial suppport. N.E.L.H. acknowledge the Research project Gaspro, financed by the research council of Norway ( 267916 ) the European Union's Horizon 2020 research and innovation programme (No 764697 ). This work also benefitted from computer resources made available through Norwegian NOTUR program, under award NN9405K.
Appendix A. Boundary treatment with immersed boundary (IB) method in foam-extend ( Jasak et al., 2014 ) In the IB method, the mesh is categorized into three types of cells called IB cells, Fluid cells or solid cells, which is shown in the Fig. A.1 a.   Fig. A1. (a) Different cells around an Immersed boundary(IB), IB cell normals, (b) Extended stencil around an IB and local co-ordinate system for Neumann boundary conditions. adopted from Jasak et al. (2014) with the permission from the authors. Velocity (Dirichlet boundary condition) of an immersed boundary cell( φ p ) is calculated using quadratic interpolation as φ p = φ ib + C 0 (x P − x ib ) + C 1 (y P − y ib ) + C 2 (x P − x ib )(y P − y ib ) + C 3 (x P − x ib ) 2 + C 4 (y P − y ib ) 2 , (A.1) and pressure (Neumann boundary condition) of an immersed boundary cell is calculated as where the coefficients C 0 , C 1 , C 2 , C 3 and C 4 are calculated using weighted least squares method in the extended stencil shown in