Kinematics and dynamics of suspended gasifying particle

The effect of gasification on the dynamics and kinematics of immersed spherical and non-spherical solid particles have been investigated using the three-dimensional lattice Boltzmann method. The gasification was performed by applying mass injection on particle surface for three cases: flow passing by a fixed sphere, rotating ellipsoid in simple shear flow, and a settling single sphere in a rectangular domain. In addition, we have compared the accuracy of employing two different fluid–solid interaction methods for the particle boundary. The validity of the gasification model was studied by comparing computed the mass flux from the simulation and the calculated value on the surface of the particle. The result was used to select a suitable boundary method in the simulations combined with gasification. Moreover, the reduction effect of the ejected mass flux on the drag coefficient of the fixed sphere have been validated against previous studies. In the case of rotating ellipsoid in simple shear flow with mass injection, a decrease on the rate of rotation was observed. The terminal (maximum) velocity of the settling sphere was increased by increasing the ratio of radial flux from the particle boundary.


Introduction
Solid particle gasification is part of many industrial processes such as combustion of carbonaceous solid fuels in a reactor or gasification of gas hydrate particles in an air lift pump [1,2]. Optimizing these processes to decrease the cost and environmental impacts requires having a profound insight about gasifying suspension behavior in flow.
The size and shape of particles are within physical properties of the solid fuels, which strongly affect observed phenomena in combustors and gasifiers. For instance, the rate of gas-solid reactions depends on both available surface area and particle size [3]. Therefore, the study of each individual gasifying particle with different size and shape helps to understand certain properties of the particle flow in large scale.
The motion of spherical and non-spherical particles has already been studied in different flow cases such as the settling of single spheres under gravity [4][5][6][7][8] and the rotation of an ellipsoid in simple shear flow [9][10][11][12][13][14]. The influence of mass injection from the surface of the solid particle has also been studied by previous authors [15][16][17] by implementing outflow on the boundary of a fixed sphere. In later works it is shown that there is a reduction in the total drag coefficient by injecting radial flux from the surface of the sphere which is correlate with mass injection ratios. However, to the best of our knowledge, the effect of mass injection on non-spherical or suspended solid particle behavior has so far remained untouched.
There are several methods that have been applied to simulate suspension particles, e.g., the finite element method [18,19] or the immersed boundary method, in original form or in combination with other methods [20][21][22]. As a requirement, all these methods must provide a proper solution to satisfy the no-slip condition on the surface of the solid particle [23]. During the last decades, the lattice Boltzmann method (LBM) has attracted much attention, and due to its ease of implementation on parallel architectures and computational efficiency, it has been applied widely to simulate hydrodynamic phenomena. Furthermore, this method benefits from a coupled timescale for both fluid and particle. Among the applications of the LBM, simulation of the suspended particles are very popular. In the first approach to this problem, Ladd [24,25] introduced the boundary condition of the moving particle by a bounce-back method on particle distribution functions. In this method, the entire domain, including the particle interior, occupied by fluid and particle inertia is influenced by the internal fluid inertia. This helps to cancel out the impulse forces exerted on the boundary due to the particle motion. A few years later, Aidun et al. [26] introduced another method in which the internal fluid inside the particle was excluded. This approach requires the so-called impulse forces due to the particle motion to remove the fluctuations in the hydrodynamic forces [26][27][28]. Furthermore, to improve the curved fluid-solid boundary in the LBM, several approaches have been suggested that utilize inter-extrapolation of the velocity or particle distribution functions [29][30][31].
In the present study, we consider the effect of radial mass flux on the motion and dynamics of a solid particle. For that reason, the rotational dynamics of the spheroid in simple shear flow and the terminal velocity of the settling single sphere are compared for the cases with and without source term. Due to the isothermal assumption of our simulation, the mass transfer is uncorrelated to the heat transfer. Therefore, an outflow condition is applied homogeneously at the solid particle surface. One can certainly argue about these assumptions. However, the aim of this work is a qualitative study of the dynamics and kinematics of the immersed particle in the presence of gasification. To avoid artifacts in our numerical results, various simulations with different domain and grid sizes have been conducted. In Sect. 4, the results of two different approaches for fluid-solid boundary are compared and validated based on previous works. Also, an accuracy study is performed for each boundary method in combination with mass injection in Sect. 4. This study helps us to identify the most suitable method for future works. The effect of ejecting mass flux on the kinematics and dynamics of the suspended particle is discussed in Sect. 5.

Lattice Boltzmann method
The LBM is an alternative method to solving the Navier-Stokes equation for flow simulation. In this method, the domain is discretized into a set of lattices. The so-called pseudo-particles are defined to be either in an equilibrium state or move through these lattices. Depending on the dimension of the domain and the number of existing probable directions for pseudo-particles to move, one can use different methods of discretization. For the current study, we used D3Q19 in which D3 stands for the dimension and Q19 shows the number of particle distribution functions for each lattice. Particle distribution functions (PDFs) denoted by f i (x, t) specify the occurrence of moving each pseudo-particle in the ith lattice direction with the velocity c i at position x and time t. The discrete velocities, c i , are defined as One method to update PDFs during time evolution is the single relaxation time based on BGK [32] approximation. In this method, the collision and streaming sub-steps are defined as in which τ and t denote the dimensionless relaxation time and computation time step, respectively. Furthermore, f i (x, t) is called the post-collision distribution function and f is the equilibrium distribution function in the D3Q19 method calculated by [33], f eq i = 3 Simulation of suspended particle

The equations of motion
The solid particle motion can be determined by translation of center of the mass, x c , and the rotation around the principal axis. The time evolution of the translational velocity, U, is obtained by solving Newton's equation in which m denotes the mass of the particle and F is the force exerted on it. On the other hand, the time evolution of the orientation of the particle is carried out by solving the Euler equation where T and Ω stand for torque and angular velocity, respectively. Moreover, I is the inertial tensor that is diagonal in the particle frame and the principal moments of inertia, the diagonal elements, are constant. Therefore, it is more convenient to solve the Euler equation in the particle frame. To transfer tensors and vectors from the particle frame (body-fixed) to the world frame (space-fixed), a so-called rotational matrix R is required (see Fig. 1). To avoid the problem of divergence in the orientational equations of motion, an alternative way is to define the rotational matrix R by quaternion parameters, instead of the Euler angles φ, θ , and ψ [34]. A parametrized quaternion is a set of four scalar quantities, q = (q 0 , q 1 , q 2 , q 3 ), with the following constraint which ensures the orthogonality of the rotation matrix: where q i s are determined by the Euler angles as Therefore, the rotation matrix based on quaternions is defined as and converts the vector v between frames, Here, the superscripts "s" and "b" denote space-fixed and body-fixed framework correspondingly. Using quaternions, the equations of motion for angular velocities can be written as ⎡ For the present study, the recent predictor-corrector direct multiplication algorithm presented by Zhao et al. [35] is applied to solve equation of motion (9) (in the body-fixed frame) and (14).

Treatments of a moving solid boundary
The velocity boundary condition at the suspended particle wall is handled by distribution functions f i . After the collision step, the distribution functions coming from solid nodes to fluid nodes are unknown. To fulfill this requirement and complete the consecutive streaming step, a bounce-back scheme [36,37] based on exchanging the momentum between solid and fluid nodes can be applied. In the bounce-back method, it is assumed that all pseudo-particles which are passed the interface and entered the solid particle, are bounced back to the fluid nodes. Therefore, for still boundaries, the values of unknown distribution functions are obtained by f¯i (x f , t) = f i (x f , t) and for moving walls by in whichī denotes the direction opposite to i and u w is the boundary velocity at x w (see Figs. 2a, b). The last term in (15) defines the amount of momentum that is added to the bounced-back pseudo-particle. In the model proposed by Ladd [24], the suspended particle boundary is presented as a stair-case zigzag boundary, passing through the midpoint node of the solid and fluid lattices (Fig. 2a). Therefore, the boundary velocity u w is calculated at the midpoint, where x c is the center of mass. This means that the non-smooth boundary requires enough resolution to avoid losing accuracy. One approach to overcome this problem is applying an inter-and extrapolation of the velocity at the boundary to determine a better approximation for the midpoint velocity [31]. Another approach to avoid the geometric discontinuity of the zigzag model is interpolation of distribution functions [30,38] at a node captured on the exact position of the particle boundary (Fig. 2b). Thus, Eq. (16) for the wall velocity changes to where α shows the fraction of intersected link between the solid and fluid nodes to the real wall (x w ).

The velocity inter-and extrapolation boundary method:
In the method proposed by Yin et al. [31] the velocity which is applied to modify the streaming step is computed by inter-and extrapolation boundary method at the midpoint (x m ) between fluid (x f ) and solid (x s ) nodes (see Fig. 2a). Using α, the interpolation or extrapolation is applied to compute the velocity at x m . In this method, the parameter α is defined as and based on that there are three cases are considered. If α ≤ 1/2 and x m is located between x w and x f , a linear interpolation gives the midpoint velocity u m (x m , t), If α > 1/2 and x m is located inside the solid domain, u m is obtained by extrapolation, in which u f and u ff are the velocity of the first and second adjacent fluid nodes to the boundary. Finally, if α = 1/2, then u m = u w which implies the original Ladd method [24]. Therefore, the modified streaming step for boundary nodes is From now on, we denote this method by I /E V .

Unified interpolation bounce-back (UIBB) boundary method:
In UIBB [30,38] the distance from fluid node to the wall, α, (see Fig. 2b) is defined as The value of the bounced back distribution function is interpolated by using the neighboring fluid lattice distribution functions, in which, on the contrary to the I /E V method, the wall velocity u w is applied to define the momentum.
Note that f¯i (x f , t), which is obtained by any of these two methods, is applied to modify Eq. (2) in the poststreaming step,

Outflow from particle surface
To simulate the gasification, an outward mass flux was implemented on the surface of the particle. The outflow mass flux is determined as 2w i ρ o c 2 s c¯i · u o · n where u o and ρ o are the velocity and density of the ejecting mass. Moreover, n is the normal vector to the surface and the subscript o denotes the outflow variables. The flux ejected from the surface of the particle is carried by the bounced-back pseudo-particle which is returning to the fluid field. Therefore, to model the ejecting mass flux for the bounce-back method, Eq. (21) is modified as and Eq. (23) is turned to This means that, using this method, the normal velocity at the particle surface is added to the wall velocity. Note here that, for three-phase simulations, the density of the mass flux is the local density of the gasifying species. However, for the current simulation it is set to the density of the adjacent fluid lattice.

Updating the force and torque
The hydrodynamic force and torque exerted on the solid particle boundary is computed based on the momentum exchange (ME) method [24,26]. The momentum transfer on each boundary node is given by The hydrodynamic force on the surface of the particle is computed by the summation of the momentum transfer on each boundary node, The hydrodynamic torque can also be obtained as Due to the particle motion, in each time step, there are some nodes that are covered by the particle which induce small impulse force on it. In the meantime, there are some nodes which become uncovered and the particle should lose the same amount of the force which is obtained from the momentum of the fluid. The requirement of this force is investigated by Aidun et al. [26] for the first time. The authors introduced the socalled cover and uncover force and torque exerted on the particle. The force and torque computed for covered and uncovered particle boundary nodes are summed up to the total force and torque. Later, Chen et al. [27] showed that regardless of the method that is used to simulate moving boundary of the non-shell model, the ME requires the impulse force to satisfy Galilean invariance. In their study, the ME method is modified by the initial momentum of the net mass transfer M and determined as follows, where M is defined as In the current study, the method introduced by Chen et al. is applied to compensate the need for the impulse force.

3D flow past over stationary sphere
The implementation of two different methods for the solid particle boundary, the midpoint velocity and unified bounce-back, are validated by a flow passing a stationary sphere in a range of moderate Reynolds numbers from 10 to 50. The definition of particle Reynolds number based on the sphere diameter d is Re = U d/ν, where ν and U are kinematic viscosity and uniform velocity at the inlet, correspondingly. The uniform inlet velocity enters the domain from the left side. An outflow boundary condition is applied at the right boundary (downstream) of the domain, and to simulate the unbounded flow, the free-slip condition is imposed for the remaining boundaries (0 ≤ y ≤ H , 0 ≤ z ≤ H ). To eliminate the influence of the boundaries on the drag coefficient, we studied different ratios H/r where r is the sphere radius and H is the distance between side boundaries. As a result, we found that H/r = 20 is suitable for current study. The center of the sphere is located at z = y = H/2 and x = 5r from the left side of the domain. The drag coefficient based on the hydrodynamic force in x-direction is often expressed as in which ρ is the density of the fluid. Figure 3 shows the comparison between drag coefficients obtained by simulation with two different particle boundary conditions. The result is compared with the relation between drag coefficient and the sphere Reynolds number expressed by Abraham [39] C d = 24 9.06 2 9.06 √ Re The black dots in Fig. 3 represent the drag coefficients reported by Nikrityuk and Meyer [40] using the immersed boundary (IB) method. The results from both methods, UIBB and I /E V , show almost the same accuracy except for particle Reynolds number 10. Furthermore, the maximum deviation from the value of C d is 12% that is obtained by relation (33) and is 8% from IB solver [40] at Re = 10. However, this error drops about 3% by increasing the resolution of the particle from 5 to 10 lattices per radius. The same effect of the particle resolution on drag coefficients is reported by Hölzer et al. [41].   . 4) with initial orientation of (φ, θ, ψ) = (π/2, π/2, 0). The major semi-axis is set to b = 12 with the aspect ratio of a p = 2. The center of the ellipsoid is placed at the center of the box with size of (N x , N y , N z ) = (96, 96, 96) (N y denotes wall distance). In this simulation, the shear rate, G = 2U/N y , and the particle Reynolds number, Re = 4Gb 2 /ν, are set to be 1/34,560 and 0.05, correspondingly. The validation is performed against Jeffery's analytical solutions [42] (Re p = 0) for angular rotation and angular rate of rotation, where c is the minor semi-axis. Figure 5 shows a good agreement with the Jeffery's solutions for both methods, I /E V and UIBB.

Settling single sphere under gravity
A settling single sphere under gravity is simulated to compare the application of I /E V and UIBB methods in presence of external force. The result is validated against the experimental data of Ten Cate et al. [7] for the particle Reynolds number (Re = u max d p /ν) 1.5 and 31.9. In the experiment conducted by [7], a spherical particle with diameter 15 mm was hanging at distance 120 mm from the bottom of a close container. The container depth and width were chosen as 100 mm, and the height was assumed to be 160 mm. To perform numerical simulation, the length and time scaling factors were obtained by so-called hydrodynamic radius, r h , and the kinematic viscosity (see [7], Table II). For the current simulation also, the domain is assumed to be a container. However, in order to satisfy mass conservation for the cases that we discuss later in chapter 5.3, the top side was prescribed as outflow boundary condition. Taking a closed container, the motion of the settling particle can be divided into three stages: acceleration, steady fall and deceleration. The steady state occurs when the drag force is in balance with the gravitational force and the particle reaches its terminal velocity. This step is observed for particle Reynolds number 31.9 within a very short time. As the particle approaches the bottom of the container, it starts to decelerate since the fluid squeezed between the particle and the bottom wall causes a force in opposite direction of the gravitational force. The driving force exerted on the sphere during the sedimentation is obtained by balancing drag force, gravity, and buoyancy, Figure 6a shows the validation of the current simulation against experimental data from Ten Cate et al. [7]. The result shows good agreement for Re = 1.5, although there is a deviation of about 12% for Reynolds number 31.9 with the resolution of 4 lattices in radius. The deviation decrease to 7.8% for the resolution of 8 lattices.
Furthermore, Fig. 6b shows that the same maximum velocity is obtained using I /E V and UIBB methods. However, I /E V gives a slightly more oscillatory result.
The results from simulation of settling sphere for Re = 31.9 reveals that final velocity is influenced by the particle resolution. In the study by Yu et al. [43], an adaptive mesh refinement method is applied for the grid resolution near the solid surface in which for the finest mesh, 60 grid points were taken across the particle diameter. Figure 7 shows a decreasing trend of the error by increasing the particle resolution, varying the radius from 4 to 32 lattice units. Note that to achieve the desired resolution of the particle, a refinement of the whole domain is required. Furthermore, the length and time scaling factors of the simulation with a resolution higher than 8 lattices per radius are calculated based on the particle radius (instead of r h ) and the kinematic viscosity. As is shown in Fig. 7, the error drops from 12 to 4%.  Error(%) Fig. 7 Resolution study of a settling sphere for particle Reynolds number 31.9

Error of computed mass flux
To validate the conservation of mass in the implemented model for gasification, the settling single sphere for Reynolds numbers 1.5 and 31.9 was considered. The mass flux, Flux comp , was computed on the outflow boundary at the top of the domain for different bounce-back methods on the solid particle boundary, different resolutions, and different particle Reynolds numbers. The result is compared with the calculated flux, Flux calc , and the errors are presented in Table 1. The error is computed by The mass injection flux on the surface of the sphere is calculated as Flux calc = ρu o A, where A denotes the surface area of the sphere. The outflow velocity u o is related to the maximum velocity of the settling sphere. This will be later discussed in Sect. 5.3. As can be seen from the table, the I /E V method gives a higher error compared to the UIBB method. This would be a direct result of the inter-and extrapolation of the midpoint velocity at the "midpoint" of the boundary lattice link. Note that in the presented model for mass injection, the outward velocity is summed up with the velocity at the particle boundary. This means that to reduce deviation of the I /E V method, mesh has to be refined more in order to have the midpoint node closer to the wall node. This is less evident at Re = 1.5 since the higher viscosity, i.e., smaller time step in simulation, could dampen the propagated error. Therefore, to decrease the error sufficiently, for Re = 31.9, more mesh refinement is required. Note that with a finer mesh (more lattice across radius) the results are closer to a real sphere. Also, the error related to the results using UIBB method decreases due to an increase in the resolution for both particle Reynolds numbers. One way to meet the accuracy tolerance without mesh refinement might be interpolating the outward mass flux over all outgoing population in each lattice. Due to a large error in computed mass flux, the results from simulation with I /E V are excluded in the following sections.

Effect on drag coefficient
The effect of mass injection on drag coefficient of a solid sphere was investigated in previous studies for Newtonian [15,16] and non-Newtonian [17] fluids. It is observed that in a Newtonian medium, applying a radial mass flux on the particle boundary reduces the total drag coefficient of the sphere. Chuchottaworn et al. [15] suggested a correlation for the ratio of the drag coefficients and mass injection rate, based on their numerical results, of In this correlation, Re denotes the particle Reynolds number and C d 0 is the drag coefficient of the sphere without mass injection. The variable a = u o /U is the injection ratio and relates the free stream velocity U (at the inlet) to the normal velocity u o at the surface of the sphere. This correlation is applicable for 1 ≤ Re ≤ 150 and −0.2 ≤ a ≤ 0.5. To observe the effect on the drag coefficient, the particle Reynolds numbers 10 and 50 were selected from the cases in Sect. 4.1. The mass injection effect on streamlines around the sphere is presented in Figs. 8a-d. The impact is more apparent at Re = 50 with a = 0.5 and can be observed by an increase in the wake length behind the sphere. The reducing trend of the drag coefficient against the injection ratio at particle Reynolds numbers 10 and 50 is plotted in Fig. 9. In both cases, we observe a decreasing trend in drag coefficient; however, a bigger drop occurs at higher Reynolds number. We also computed the ratio C d /C d 0 and compared it with (38) (see Fig.  9). For a = 0.1, the error is about 1% for both Reynolds numbers. However, the largest deviation for both set of data is related to the highest ejection ratio. The deviation is about 4% for Re = 10 and about 16% for Re = 50. To study domain size dependence, the simulation is repeated with a = 0.5 and for a 20% larger domain in spanwise direction (y and z). In this case, we observed that the effect was less than 1%. Therefore, the explanation of the error could be related to the nature of the lattice Boltzmann method and follow the same reasoning discussed in Sect. 4.4.

Effect on period of rotation
In this section, the effect of ejecting mass on the dynamics of the rotating particle is studied by adding the outward flux to the surface of the ellipsoid. The simulation settings are the same as the case presented in Sect. 4.2. The ejection rate is obtained by u o = aGd where d = 2b and a varies from 0 to 0.5. To satisfy  Fig. 9 The trend of the decreasing drag coefficient for the flow passing by a fixed sphere, as a function of injection ratio a mass conservation, we substitute the periodic boundary condition with the outflow boundary condition only in velocity direction (x direction). This increased the period of rotation around 1.4%. We note here that, for all the cases with mass injection, simulation with non-periodic boundary condition and a = 0 is taken as the reference.
The streamlines in Fig. 10 show the effect of mass ejection from the surface of the ellipsoid on the surrounding flow field. Figure 10 (right) shows that not only the background flow (continuing flow region) is pushed toward the walls but also, due to the outward flux with the factor of a = 0.5, the reverse flow regions in left and right are displaced. This may conclude that the mass injection causes a reduction of shear rate close to the rotating ellipsoid and increases the period of rotation. To indicate this, for the set of simulations with different ejecting factors, the periods of rotation are computed and compared together. In addition, to investigate the influence of the confinement, the simulations were repeated for the domain size of (N x , N y , N z ) = (96, 120, 96). However, to avoid the effect of time resolution, the shear rate remained unchanged. Figure 11 shows direct relation of the period of rotation and ejection ratio for both domain sizes. Despite of the average 5% difference in value, the results for different box sizes revealed qualitatively similar behaviour.
Rosen et al. [14] have studied the dynamics of prolate spheroids and discussed about different equilibrium rotational states resulted by initial orientation and increasing particle Reynolds number. We remark that the Ny=96 Ny=120 Fig. 11 The trend of increasing rotational period by increasing the injection ratio for two different wall distances N y = 96 and N y = 120 present simulation result performs a tumbling state, in which the particle symmetry axis is located in a flow gradient plane (θ = π/2). In addition, based on the discussion in [14], in a tumbling motion of the prolate ellipsoid, the particle inertia has a dominant inertial effect. In the current study, since there is no changes in the particle Reynolds number and the Stokes number (St = αRe, α is density ratio), the particle inertia does not differ with mass injection. This may support the previous discussion that the reason for an increase in the period of the rotation could be the effect of the mass injection on the flow around the rotating ellipsoid.

Effect on final velocity of settling sphere
The final velocity of a sedimenting sphere with mass injection is compared with simulation results in Sect  Table 1) is shown in Fig. 12.
The maximum velocities for different injection ratios are scaled by the terminal velocity of the case without mass injection. To understand the effect of wall friction on the final velocity for the cases with mass injection, the simulations were repeated for wider domains (by increasing distance between walls). We observed that, for both particle Reynolds numbers, u t /u 0 t > 1. This is an indicator of correlation between drag reduction and increasing the ejection rate (see Sect. 5.1). However, for the original domain size at Re = 1.5, the terminal velocity does not increase at the same pace as the one for the wider domain. This means that the mass ejection from the sphere surface may increase the wall resistance and consequently compensate the drag reduction effect.
Note that, as Fig. 13b indicates, at particle Reynolds number 31.9, the growth rate of u t /u 0 t becomes slower for a < 0.3. The maximum value of u t /u 0 t is attained at a = 0.3, and thereafter, the slope of the curve becomes negative. Beside the wall friction effect, another reason for this behavior could be the upward force in opposite direction of the gravitational force, which provokes deceleration of particle velocity. The flow that is bounced back from the bottom wall originates a resistance which gets stronger by increasing the injection rate. The effect of upward force is more apparent at higher Reynolds number where the particle advection time is smaller. Moreover, with the domain height fixed for all injection ratios, one can argue that increasing the settling velocity with increase in ejection ratio may postpone the steady fall phase. This means that the particle reaches to decelerating stage before achieving its terminal velocity. For clearer understanding, we undertook two separate simulations with 50% wider and 50% longer domains (in gravity direction) for Re = 31.9. The result in the Fig. 13b indicates that the effect of increasing wall distance is more noticeable than increasing the domain height. Therefore, we can conclude that wall friction is the reason for declining in growth rate of u t . However, there is a maximum increase of 4% in u t /u 0 t for the longer domain compared to the original domain size.
Comparing Figs. 13a and b shows a stronger increase in the final velocity at Re = 31.9. To clarify this observation and to interpret the effect of mass injection on the transient behavior of the settling sphere, we take into consideration the momentum diffusion timescale and the injection timescale. The former timescale is defined as τ ν d 2 /ν, and the latter as τ inj d/u o , in which d stands for the particle diameter and ν is the kinematic viscosity. The momentum diffusion timescale estimates the time for diffusion of momentum into the fluid over a distance of one particle diameter. Based on this definition, the diffusion timescale at Re = 1.5 is τ ν = 0.59 and at Re = 31.9 is τ ν = 3.72. Table 2 presents mass injection timescales based on different injection ratios. The comparison between diffusion timescale and the data set in Table 2 indicates that the ratio τ ν /τ inj at Re = 1.5 is smaller than one; however, this ratio for Re = 31.9 is greater than one. This means that at Re = 1.5 there is time for the injected mass to be adjusted into the surrounding fluid. In contrast at Re = 31.9, the momentum diffuses one diameter into the fluid at slower rate than the injection, which means that the surrounding fluid is disturbed stronger.
In Fig. 14, the effect of ejecting gas on the fluid flow around the settling sphere is shown. This change is apparent in both the front and rear parts of the sphere. In addition, it can be seen that due to the mass injection the circulation of the fluid is pushed away from alongside of the sphere.

Conclusion
The lattice Boltzmann method was used to simulate rigid immersed particles by applying two methods for particle boundary condition, the velocity inter-and extrapolation boundary method and unified interpolation bounce-back. The results of 3D flow cases, flow past by a fixed sphere, a rotating ellipsoid in a simple shear flow and the sedimentation of a single sphere under gravity, were validated against previous studies. The same accuracy was observed for both the UIBB and I /E V methods. A model was presented to simulate gasification by implementing mass injection from the surface of the particle. The accuracy of the model was investigated by comparing the calculated mass flux against the one from the numerical simulation results. For that purpose, the error between calculated mass flux on the surface of the sphere and obtained mass flux from the simulation of the settling sphere (gathered on the top of the domain) was studied. The highest deviation was observed in the result of the simulations with the I /E V method. Therefore, this method was not found suitable to simulate the gasification with our suggested model. A dependency on the resolution was observed in the error. This can be a direct result of simulating the sphere using a Cartesian equidistant mesh in lattice Boltzmann method. To investigate the effect of implementing gasification on the kinematics and dynamics of the particle, the mentioned flow cases were simulated under the condition of ejecting mass from the surface of the particle. The reduction of the drag coefficient of the sphere due to mass injection was validated against data from suggested correlation by Chuchottawornin et al. [15]. The results show a deviation for higher injection ratios and particle Reynolds number 50.
Moreover, the influence of the outflow on the dynamics of the suspended particle is investigated. It is shown that the source of the flux on the boundary has an effect on both rotational and translational motions of the particle. The effect of mass injection from the rotating ellipsoid surface in simple shear flow was detected by the reduction of the rotation speed. Therefore, we observed an increase in the period of the rotation compared to the case with ejection ratio 0. Also, ejecting mass flux caused an increase in the final (maximum) velocity of the settling sphere in a 3D rectangular box. However, for the case with higher Reynolds number, a smaller increase in the maximum velocity was observed.
To conclude, we claim that the current model is able to simulate suspended spherical and non-spherical particles with gasification. However, some improvements might be required to increase the accuracy of the method particularly for higher Reynolds numbers and higher injection ratios.