CHEMICAL ENGINEERING TRANSACTIONS VOL. 70, 2018 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Timothy G. Walmsley, Petar S. Varbanov, Rongxin Su, Jiří J. Klemeš Copyright © 2018, AIDIC Servizi S.r.l. ISBN 978-88-95608-67-9; ISSN 2283-9216 Simulation of Gas-Solid Flow in Quasi-Two-Dimensional Fluidized Bed by Immersed Boundary Method Lívia Gyurik*, Attila Egedy, Zsolt Ulbert University of Pannonia, Dep. Process Engineering, Egyetem St. 10, 8200 Veszprém, Hungary gyurikl@fmt.uni-pannon.hu Gas-solid two-phase flow occurs in multiple fields of industry. Therefore, a study researching the phenomena of this type of systems is highly adequate. In the present work subsonic, compressible flow equation in conservative Euler form is solved numerically on PC by a two-step finite difference method, which is validated by the analytical solution of Sod’s shock tube problem. Fluid-solid interaction is described by Immersed Boundary Method, where body force is calculated by direct forcing algorithm. Two-dimensional simulation of a single solid circle shaped particle motion in the gas flow is calculated and compared to the trajectory of a particle recorded by a high-speed camera measurement system in a laboratory-scale quasi-two-dimensional fluidized bed. Image processing algorithms are implemented to track the motion of a chosen single. Comparison of experimental measurement and simulation resulted in the same magnitude of particle velocities. The possible reasons of differences and opportunities of improving the results are discussed in the study. The results are promising, and the presented preliminary simulator can be a useful tool in supporting industrial design and operation after further developments. 1. Introduction Gas-solid two-phase flow occurs in many processes in chemical and pharmaceutical industries e.g. granulation, mixing, pneumatic transport and fluidization. Fluidized beds are widely used in processes like fluid catalytic cracking (John et al., 2017), particle separation (Weber et al., 2017) and drying (Idakiev et al., 2017). In gas-solid flows, besides the solid phase, gas phase also plays an important role, and the efficiency of the processes is determined by the gas-solid interaction. Multiphase interactions have huge energy demand; however, costs can be reduced by optimization. The measurements in gas-solid multiphase flow systems are limited because they can influence the flow behavior of the system. In this way, simulations provide a good possibility to understand the behavior of the gas-solid flow system and help to determine the optimal operating parameters. Optimal design of equipment can provide energy, time or size gain (Jurena et al., 2017). Implementation of changes in the geometry of experimental setup is often limited due to its cost and material needs. The intention for environmental protection also facilitates the research of gas-solid two-phase flows. Brachi et al. (2017) studied biomass to explore potential feedstock for bioenergy production by a deeper understanding of segregation and fluidization behavior. These points emphasize the benefits and relevance of modeling and simulation. To ensure the applicability of the models and simulations, experimental validation should be performed. There are different approaches to mathematical modeling of gas-solid two-phase flows. The first approach is the Euler-Euler method and application of the Two-Fluid models (TFM) where the solid phase treated as fluid as well as the gas phase. In this case, results are less detailed, but calculations can be performed with less computation cost. Computational Fluid Dynamic (CFD) methods and software can be applied to calculate Euler-Euler models. The second approach is called Discrete Element Method (DEM) which treats the solid phase as individual particles. It computes their collision with each other and the wall, and the movement is calculated by the sum of contact forces emerge on the particle. This approach gives detailed information about the solid phase, but it needs large computational capacity. The method was used to simulate many applications, for instance, examining solid phase mixing in a fluidized bed (Oke et al., 2016), modeling motion of pharmaceutical tablet shapes in a film coating pan (Ketterhagen, 2011) or simulation of spouting of corn DOI: 10.3303/CET1870135 Please cite this article as: Gyurik L., Egedy A., Ulbert Z., 2018, Simulation of gas-solid flow in quasi-two-dimensional fluidized bed by immersed boundary method , Chemical Engineering Transactions, 70, 805-810 DOI:10.3303/CET1870135 805 (Ren et al., 2012). Simulation of two-phase flow and calculation of the motion of solid particles by DEM can also be performed by CFD software packages using user-defined functions. Direct numerical simulation, for example, Immersed Boundary Method (IBM), is a further way to model solid-fluid interaction. The main advantage of direct numerical simulations that they are capable the most detailed description of gas-solid flows. In the Immersed Boundary method, first introduced for describing blood flow in the heart where valves acted as a moving boundary (Peskin, 1972), structured calculation mesh can be applied, and the inner boundaries (e.g., particles) can be virtually defined using body forces and extra force term added to the momentum equation. Structured (or uniform, Cartesian) mesh is often called Eulerian grid while the node points of the immersed boundary referred as Lagrangian points. These two coordinate systems are independent, but an interpolating function gives the relation between locations and flow variables. In our work, the fluid phase is considered as a real compressible gas phase. Therefore compressibility brings some additional complexity into the calculation. Flow problems can be described mathematically by partial differential equations. The main methods of solving these equations numerically are the finite element, the finite volume, and the finite difference methods. One of the possible ways of solving the conservation form of fully compressible flow equations of the gas phase is by using the two-step MacCormack’s finite difference method (MacCormack, 1969). The method employs a predictor and a corrector step which give second-order accuracy both in time and space. A physical flow field with high gradients and discontinuities, such as shock waves, does not exhibit oscillations in the vicinity of discontinuities. On the other hand, the presence of discontinuities makes the calculation of the gas flow a difficult problem to solve because the conservative second-order numerical methods generate strong oscillation around the discontinuous solution which can cause errors in the simulated gas flow field and fluid drag force calculation. Thus, the numerical approach must effectively handle steep spatial fronts. To handle the steep spatial fronts, it is natural to apply modern shock-capturing numerical methods, as for instance Total Variation Diminishing (TVD) schemes introduced by Harten (1983), for the convective part of the conservation laws. Applying a TVD scheme, we used the explicit TVD-MacCormack predictor-corrector method in which a second corrector step, addition of an appropriate conservative dissipation term is introduced (Yee, 1989). Sod’s shock tube problem serves commonly as a validation of the applied numeric method (Sod, 1978). Its analytical solution is known, so the results of the applied solver can be compared to the exact solution. Out of those processes where gas-solid two-phase flow occurs, we focus on fluidization process as a case study in this work. In fluidization systems, the injected gas stream keeps the particles in a fluidized state. In our simulation study, we simulate the particle motion in a quasi-two-dimensional laboratory-scale fluidized bed. Our model is based on the compressible flow equation, and IBM modeling strategy for the fluid-solid interaction. A code is developed in MATLAB environment to solve the flow equations. The results are experimentally validated with the use of a high-speed video camera measurement system to capture the real motion of the particles and be able to catch particle trajectories by image processing and compare to the simulated particle motion. 2. Materials and methods 2.1 Flow equations As the gas phase possesses the ability of flowing, it is called fluid as well as liquids. The difference is in their compressibility features, i.e., liquids are incompressible, while gas phase is compressible. Viscosity gives the other difference since in case of gas flow viscosity is usually omitted, while it plays an important role in case of liquid flow. These differences appear in the fluid equations. For compressible, viscous flow usually the Navier-Stokes equation is used, while for inviscid flow the Euler equation is applied. These equations describe the change of flow variables in time and space. The present research uses the conservative form of Euler equations as it is written in Eq(1) - (3) for one-dimensional case. + ( ) = 0 (1) ( ) + ( ) + = 0 (2) + + + + ( ) = 0 (3) Where ρ is density, v is velocity, p is pressure, e is internal energy, t is time, x is space. 806 Numerical solution of these types of partial differential equations can be performed by several numerical methods. The finite difference numerical methods, for example, Lax-Wendroff (Lax and Wendroff, 1960) and MacCormack’s technique (MacCormack, 1969), are frequently used methods in solving hyperbolic equations. Both methods are second-order accurate in time and space, but the implementation of the latter is easier, as it uses only first-order differentials. Eq(1) – (3) are solved by MacCormack’s two-step predictor-corrector method. The inlet boundary condition was defined by a constant air inflow velocity 3.33 m/s in accordance with the experimental equipment, while constant pressure boundary is applied at the outlet. At the vertical walls, the no-slip condition is used. 2.2 Validation of the numerical method Flow problems with analytical solutions can serve as validation cases. One of the often used test problems is the so-called shock tube problem. In Sod’s example (Sod, 1978), a diaphragm in the middle of a tube separates the two areas with different pressure and density conditions. After removing the diaphragm, a flow spreads in both directions. In the case of general space-centered numerical methods, unexpected oscillations emerged beside the steep spatial fronts (Figure 1a). These are not observable in physical flow; therefore it implies a numerical problem. There are a lot of modern shock-capturing numerical methods to handle steep spatial fronts, one of them are Total Variation Diminishing (TVD) schemes. It attenuates the oscillations effectively as Figure 1b shows. Figure 1: The effect of TVD on the numeric solution. a) solution without TVD, b) solution with TVD 2.3 Fluid-solid interaction Gas-solid two-phase flow is simulated by Immersed Boundary Method. This method requires the introduction of an extra force into the flow equation to modify the fluid flow in accordance with the solid body surface (Eq(4)). ( ) + ( ) + + = 0 (4) Where f is the extra force. The applied forcing scheme that is used to calculate the body force follows the direct forcing method as described in (Fadlun et al., 2000). The main idea of this method is to set the value and direction of f to balance momentum at those grid points where the solid immersed object takes place. Coupling between Eulerian and Lagrangian grids is carried out by an interpolation stencil. The effect of a body force point on a fluid grid point is weighted by the interpolation stencil. On the other hand, the effects of fluid properties of fluid grid points on a solid body point are also weighted by the same stencil using weighting function given by Eq(5). ( ) = | | , 0 | | 0, | | (5) Where w is weight, d is the distance between body surface point and fluid grid point, and h is the grid size. Position update of the particle is calculated by force emerging along the immersed boundary. Pressure values around the object can be obtained from the flow equations, which are divided by the area, and multiplied by the normal vector of the given area. Acceleration, velocity, and position are calculated from the force by numerical integration, respectively. Calculations are performed on a Dell Optiplex 790 PC with 16 GB memory. 807 2.4 Experimental setup and video processing A laboratory-scale fluidized bed is built up for experimental validation. The quasi-two-dimensional equipment (Figure 2) is 14.5 cm wide, 80 cm high, and 1.5 cm deep. The speed of air inflow is adjustable and measurable, and pressure values at the left side wall at different heights of the equipment can be logged by the use of Advantech ADAM-5000L/TCP data acquisition and control system. The fluidized bed is filled with approximately 11,000 perfect sphere shaped metal particles with a size of 3.5 mm in diameter and weight of 0.067 g. The video for validation is captured by Marathon Ultra-high-speed camera measurement system including Optronics CL600x2 camera with 500 frames per seconds providing raw data for image processing. In the image analysis, a suitable region of interest is cut off from the sequence of frames where an upward motion of a single sphere is observable. Particles in the greyscale picture are detected by a shape recognizing algorithm. Figure 2: The laboratory-scale quasi-two-dimensional fluidization equipment 3. Results and discussion Numerical solution of the flow, Eq(1) – (3), is performed by the second-order MacCormack finite difference method which is validated by comparison to the analytic solution of Sod’s shock tube problem. According to (Sod, 1978) initial condition for density is 1 kg/m3 in the first half of the tube and 0.125 kg/m3 in the second half. After removing the diaphragm in the first time step, a sharp front starts to spread. In the numerical solution, erroneous oscillation occurs at steep spatial fronts which not follow the real physical variation of variables (Figure 1a). To handle this numerical problem, Total Variation Diminishing scheme is used to smooth oscillations. The method successfully eliminates the numerical oscillations which can be noticed in Figure 1b. In order to simulate two-dimensional particle motion in air flow, the Immersed Boundary Method was used. The compressible air flow is calculated by solving the conservation form of two dimensional fully compressible Euler equations by the validated MacCormack’s numerical method. The solid particle is represented in the calculation domain by its Lagrangian grid. The calculation of the effects of a solid body on the gas flow is carried out by adding an extra force term to momentum equations as shown in Eq(4). The effect of a solid body on the gas flow is shown in Figure 3. The gas flow around the particle body is indicated by the velocity vectors. The grid size of gas flow calculation was chosen to 0.001 m, 0.0008 m, 0.0006 m and 0.0004 m respectively, while time step was calculated according to Courant-Friedrichs-Lewy condition (Courant et al., 1928). Grid size advantageously influences the accuracy of the solution and the computational time. Refining the calculation mesh increases the computational time and decreases the velocity of the particle at simulation time 0.014 s respectively (Table 1). Using a captured frame taken by the high-speed camera, we selected a particle with upward motion in its region of interest. The selected particle was tracked, and its trajectory was compared to the motion of the simulated particle. In the region of interest of the captured frames, the circle shapes are recognized due to an image processing algorithm, and the position of the selected particle (dark filled particle in the white colored rectangle area, Figure 4) is logged. In the simulation, 14,083 time steps were calculated to reach the simulation time 0.008 s, and 24,905 time steps to the simulation time 0.014 s. Pictures are taken in every 808 0.002 s by the high-speed camera. The position of simulated and captured particle at three different times is shown in Figure 4. The chosen grid size is 0.0004 m, the finest out of the examined versions. The velocities of the particle obtained by the simulation are higher compared to the results of the captured video frame (Table 2). The reason of the difference is that only a single particle is calculated in the simulation, and the inlet gas velocity of calculation domain is overestimated by the calculated experimental superficial gas velocity. The setting of the boundary condition of inflow also contains inaccuracy since determining the velocity at a region of interest in the experimental equipment is inhibited. In reality, pressure is not reflected on the outflow boundary. However, the model works with reflective, constant pressure boundary condition. Adding a non-reflective scheme will enhance the realism of the simulator. Oscillations at steep spatial fronts still exist as TVD has not built into the two-dimensional version yet. Figure 3: Immersed boundary in the flow field. Grid sizes: a) 0.001 m, b) 0.0008 m, c) 0.0006 m, d) 0.0004 m Table 1: Mesh dependence of the simulation Grid size Computational time Time steps Elapsed time Total change in position Velocity 0.001 m 18.8 min 10,000 0.0140 s 0.024 m 1.72 m/s 0.0008 m 30.5 min 12,500 0.0142 s 0.020 m 1.40 m/s 0.0006 m 61.3 min 16,500 0.0141 s 0.015 m 1.05 m/s 0.0004 m 260.4 min 25,000 0.0141 s 0.010 m 0.69 m/s Figure 4: Three snapshots of the simulated, and the video processed particle (dark filled particle in the white colored rectangle area) motion in a region of interest at 0 s, 0.008 s, and 0.014 s. Table 2: Comparison of changes in position and velocities in case of simulation (S) and experiment (E) Time interval Change in position(S) Change in position(E) Calculated vel.(S) Calculated vel.(E) Error 0-0.008 s 0.0056 m 0.0039 m 0.7 m/s 0.49 m/s 41.7% 0.008-0.014 s 0.0041 m 0.0029 m 0.68 m/s 0.48 m/s 42.9% 0-0.014 s 0.0097 m 0.0068 m 0.69 m/s 0.49 m/s 42.2% 809 4. Conclusions In this study simulation of gas-solid flow in quasi-two-dimensional fluidized bed was performed. Flow equations were numerically solved by the two-step MacCormack’s method, which is second-order accurate both in time and space. The solution produces unreal oscillations in the vicinity of discontinuities as all second or higher order space-centered finite difference methods, which make the calculation of the gas flow a difficult problem to solve. A TVD scheme was applied successfully to eliminate numerical oscillations in a one- dimensional problem. Implementation of TVD scheme into the two-dimensional case will be essential to enhance the reliability of the simulator. Compression wave reflection due to constant pressure boundary condition at the outlet in case of compressible flow obstructs the formation of the stationary flow field. By applying a so-called non-reflective boundary condition, more realistic flow pattern can be formed especially in case of compressible flows. Therefore including a non-reflective boundary condition in the present model is important. The interaction between gas flow and the solid particle was handled by Immersed Boundary Method with direct forcing scheme. The simulations results show that particle velocity determined by the image processing and the simulation is in the same order of magnitude. However, refinements are definitely needed to improve the accuracy of the model. Forcing schemes can be refined by the use of a smoother interpolation stencil, and other approaches to calculating the body force. This will also help to determine position update of solid particle more accurately. Obviously, the simulation will be extended to multiple particles. Acknowledgments The research was supported by EFOP-3.6.1-16-2016-00015 Smart Specialization Strategy (S3) - Comprehensive Institutional Development Program at the University of Pannonia to Promote Sensible Individual Education and Career Choices project. References Brachi P., Chirone R., Miccio F., Miccio M., Ruoppolo G., 2017, Segregation and fluidization behavior of poly- disperse mixtures of biomass and inert particles, Chemical Engineering Transactions, 57, 811–816. Courant R., Friedrichs K., Lewy H., 1928, Über die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen, 100 (1), 32–74. Fadlun E.A., Verzicco R., Orlandi P., Mohd-Yusof J., 2000, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, Journal of Computational Physics, 161 (1), 35– 60. Harten A, 1983, High Resolution Schemes for Hyperbolic Conservation Laws, Journal of Computational Physics, 49 (3), 357–393. Idakiev V.V., Lazarova P. V., Bück A., Tsotsas E., Mörl L., 2017, Inductive heating of fluidized beds: drying of particulate solids, Powder Technology, 306, 26–33. John Y.M., Patel R., Mujtaba I.M., 2017, Modelling and simulation of an industrial riser in fluid catalytic cracking process, Computers & Chemical Engineering, 106,730–743. Jurena T., Dohnal M., Hajek J., 2017, Simulation of multiphase flow in stirred tank of nonstandard design, Chemical Engineering Transactions, 61, 625–630. Ketterhagen W.R., 2011, Modeling the motion and orientation of various pharmaceutical tablet shapes in a film coating pan using DEM, International Journal of Pharmaceutics, 409 (1–2), 137–149. Lax P., Wendroff B., 1960, Systems of conservation laws, Communications on Pure and Applied Mathematics, 13 (2), 217–237. MacCormack R., 1969, The effect of viscosity in hypervelocity impact cratering, American Institute of Aeronautics and Astronautics, 69-354. Oke O., Van Wachem B., Mazzei L., 2016, Lateral solid mixing in gas-fluidized beds: CFD and DEM studies, Chemical Engineering Research and Design, 114,148–161. Peskin C.S., 1972, Flow patterns around heart valves: a numerical method, Journal of Computational Physics, 10 (2), 252–271. Ren B., Zhong W., Chen Y., Chen X., Jin B., Yuan Z., Lu Y., 2012, CFD-DEM simulation of spouting of corn- shaped particles, Particuology, 10 (5), 562–572. Sod G.A., 1978, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, Journal of Computational Physics, 27 (1), 1–31. Weber J.M., Stehle R.C., Breault R.W., De Wilde J., 2017, Experimental study of the application of rotating fluidized beds to particle separation, Powder Technology, 316,123–130. Yee H.C., 1989, A class of high resolution explicit and implicit shock-capturing methods, NASA TM 101088. 810