Numerical modelling of flow field within a packed bed of granular material

Computational domain discretization is the pre-processing stage in CFD analysis and it highly contributes to reducing the level of solution error as well as increasing the numerical stability of a model. The discretization of domain representing a packed bed of granular material is a demanding task due to the topology of the model consisting of spherical granules contacting tangentially and narrow regions in the direct vicinity of the contact point. Therefore the aim of the carried out research is to define guidelines to effectively discretize the computational domain representing the packed bed of granular material and validate it using Particle Image Velocimetry (PIV). Two methods of contact point representation between granules (contact point extension, granule radius reduction) and two mesh types (tetrahedral, polyhedral) are examined within the paper. Differences between the analyzed cases are discussed and compared to experimental research results. The obtained results prove that the applied method significantly affects the calculated porosity but does not affect the flow field considerably. Polyhedral mesh turned out to be preferable as the numerical model converges within substantially lower number of iterations and the obtained results are not as affected by the numerical diffusion as in the case of tetrahedral mesh.


Introduction
Granular materials are commonly used in many industries and applications; therefore they are the focus of attention of numerous research efforts carried out utilizing both experimental and numerical techniques [1][2][3][4][5][6][7]. One of the latter is Computational Fluid Dynamics (CFD), which has proven its applicability as an efficient research tool for numerous industrial sectors and can also be successfully applied in flow analysis within packed beds of granular material [8,9].
The initial and crucial step in CFD analysis is the discretization of the computational domain (meshing) as it influences not only the result accuracy but also the calculation time and numerical stability of the model [10]. Discretization consists of dividing the analyzed geometry into numerous control volumes commonly called elements or cells.
Three mesh types depicted in Figure 1 are used in CFD applications: hexahedral (HEX), tetrahedral (TETRA) and polyhedral (POLY). The first (HEX)is characterized with low numerical diffusivity, but it is difficult or even impossible to be generated in the case of complex geometries. The second (TETRA) is relatively easy to be generated, but it introduces significant errors resulting from its high diffusivity. Therefore, the POLY mesh is a beneficial compromise as it combines acceptable numerical diffusivity with ease of generation. The benefit of POLY results from the fact that each individual cell neighbors with many others. In contrast, TETRA neighbors with only four elements. POLY contributes to better gradient approximations and in consequence to more reliable CFD analysis results . Moreover POLY cells are less sensitive to stretching in comparison to TETRA cells which improves the numerical stability of the model. The mesh quality is significantly important in the case of modelling a flow field in packed beds of granular materials because of the specific topology of the domain, where individual granules of spherical shape contact each other tangentially at a single point. It leads to mesh elements of extremely low quality because the cells in direct vicinity of the contact point are highly skewed. Therefore the aim of the carried out research is to define guidelines to effectively discretize the computational domain representing the packed bed of granular material.

Experimental research
The experimental research object consisted of three identical spheres of 30mm diameters. The spheres represented granules positioned inline and tangentially contacted. Each of the granules was mounted on a separate post of adjustable height attached to a wooden stand as depicted in Figure 2. The experimental research with the application of Particle Image Velocimetry (PIV) was conducted in a low-speed, closed-circuit wind tunnel. The measuring section of dimensions 0.5m x 0.5m x 1.5m was adapted for optical, non-intrusive measurement methods such as PIV (Figure 3). The free stream velocity of the tunnel ranged from 0.1 to 62 m/s and longitudinal free stream turbulence intensity was less than 0.5% as described in [11]. The wind tunnel was equipped with systems designed to control temperature and relative humidity. The utilized PIV system consisted of 5.5Mpx sCMOS camera and Nd:YAG laser emitting a pair of light pulses of energy equal to 0.2J each. Flow was seeded with ~1μm droplets of Di-Ethyl-Hexyl-Sebacate (DEHS). The seeding particles were introduced into the flow in stable concentration adjusted during the preliminary calibration and setting.
The laser light sheet illuminated plane crossing through the centers of the investigated granules. Double frame images were recorded at 15Hz frequency rate and analyzed with the adaptive correlation algorithm. Total number of double frame images taken during the measurement was 250. The research object was mounted centrally along the longitudinal axis of the wind tunnel's measuring section as depicted in Figure 3 and Figure 4. Instantaneous velocity fields were measured within the area marked by the red dotted line in Figure 3 and located in the plane coincident with the centers of the granules. The inflow velocity was equal to 1.81m/s, which corresponds to Re number related to the granule's diameter equal to 3700. The PIV research provided the velocity vector fields with two velocity components in stream-wise and span-wise directions.

Computational domain discretization and numerical model parameters
The pre-processing stage in CFD analysis of granular material domain is a challenging task due to the single-point contact between the two individual granules and very narrow space to be meshed in the direct vicinity of the contact point ( Figure 5a). It leads to extremely low-quality mesh elements and in consequence, to numerical instability of the calculations. The common approach to overcome this issue is to reduce the radius of the modelled granules (Granule Radius Reduction -GRR) while maintaining their positions ( Figure 5b). This unfortunately leads to the disturbance of flow and heat transfer in the domain. Another approach investigated in [8,[12][13][14] is to extend the contact point between the two individual granules to a cylindrical volume as shown in Figure 5c (Contact Point Extension -CPE). This approach allows generating a high-quality mesh while maintaining contact between granules. Another important issue influencing the CFD results and associated with mesh is the type of mesh elements used to discretize the domain [15]. As was mentioned earlier, HEX elements are characterized by low numerical diffusion, but the geometric complexity of the packed bed of granular material makes it impossible to apply structured HEX mesh in such a case. Therefore TETRA and POLY elements are analyzed in the paper as they can easily by applied to describe even very complicated geometries [16]. In consequence, four cases were analyzed within the numerical research: two approaches (GRR, CPE) to increase the mesh quality in the contact region and two different mesh types (TETRA, POLY).
The computational domain depicted in Figure 6 represented one quarter of the experimental research object (Figure 2) placed in a virtual wind tunnel domain of dimensions based on best practice guidelines defined in [17,18]. The geometry was created using SolidWorks with ANSYS plugin, which allows the bi-directional associativity between CAD and mesher. It was prepared in two configurations. The first represents a section of the packed bed with granules of reduced diameter (GRR method) -the reduction of granule radius was defined as (where: r -granule radius) on the basis of findings described in [9,12]. The second configuration represented the full-size granules with contact point extended to the cylindrical volume (CPE method) of radius equal to . Such radius was defined as optimal in [12].
The mesh was generated in ANSYS 18.1 with Workbench Poly Meshing for Fluent 1.8 ACT extension. The mesh sizes on the granule walls were set as with growth ratio equal to 1.05. The volume mesh size was set to . Such settings were consistent with guidelines provided in [19]. The identical computational model (except for the mesh) was used for all cases analyzed within the research. The solver was configured as pressure-based and analysis was performed for a steady state. Standard k- viscous model was applied.
Velocity-inlet boundary condition type was assigned to the inlet of the computational domain (marked green in Figure 6) with the velocity magnitude corresponding with the experimental research and equal to 1.81m/s. The outlet of the computational domain (marked blue in Figure 6) was defined as pressure outlet with gauge pressure of 0Pa (operating pressure -101325Pa), backflow turbulent intensity was set to 5% and backflow turbulent viscosity ratio was set to 10. Symmetry was assigned to the section planes indicated yellow in Figure 6.
Pressure-velocity coupling by SIMPLE algorithm was used as a solution method. This algorithm uses a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field.
The numerical calculations were performed with ANSYS Fluent 18.1 on Windows 10 workstation equipped with Intel i7 4510U CPU 2.00GHz and 16GB of RAM. The solver was operating in a single mode.

Results and Discussion
The experimental research results obtained with PIV concerning time-averaged flow field in the vicinity of the research object are depicted in Figure 7. The boundary layer close to the granules surface is clearly visible although its thickness may be misrepresented due to the resolution of the measuring system equal to 0.712mm for the analyzed region. Moreover slight differences concerning flow patterns in the two regions between the adjacent granules are noticeable.
Another PIV measurement series was performed in order to focus on the flow very close to the contact point between the adjacent granules. Limiting the examined region, resulted in improved measurement resolution equal to 0.085mm. The results concerning the vector plot, as well as time-averaged flow fields in the region between the first and the second granule are shown in Figure 8. The velocity field is presented as velocity magnitude (Velocity), stream-wise velocity (U) and spanwise velocity (V). As it can be observed, the flow velocity magnitude at the distance of 3mm (20% of granule radius) upwards from the granules' contact point does not exceed 0.02m/s which is approx. 1% of the inflow velocity.
Therefore, it can be concluded that CFD numerical analysis can be effectively performed with the CPE method of discretizing the granular material domain because even in the case of relatively large cylindrical volume radius equal to 3mm, the flow is not disturbed significantly. In contrast, the GRR method affects the main flow stream because of the reduction of the individual granule size.
In addition, the adulteration of the investigated granular material porosity is at the level of 0.03% in the case of CPE and 3.26% in the case of GRR, which directly points to the CPE method as more beneficial. As was mentioned earlier, the CFD analyses were performed for four configurations i.e. two discretization methods and two mesh types. Mesh statistics including minimum orthogonal quality (ranges from 0 to 1, where values close to 0 correspond to low quality) and maximum aspect ratio (1 is ideal aspect ratio) for all analyzed cases are listed in Table 1. The above mentioned quality indicators point to a very good quality mesh in the case of POLY and good quality mesh in the case of TETRA. The mesh quality is not affected significantly by the applied method of granule contact representation but only the POLY mesh type indicated preferable results.
The model was assumed as converged where all residuals fell below the value of 1e-3. As can be seen in the last row in Table 1, the number of iterations necessary to assure convergence was almost 2.5 times higher in the case of TETRA mesh in comparison to POLY mesh. It leads to a conclusion 6 1234567890 ''""  There are no noticeable differences in the flow field near the granules regarding the applied method of contact point representation but a decrease of granule size in the case of GRR method slightly influences the main flow above the granules. The applied mesh type seems to impact the results more significantly as lower numerical diffusion of POLY mesh is reflected in less velocity fluctuations particularly along the vertical axis of the computational domain. This issue is the most significant after the third granule. Therefore the results obtained with POLY mesh are more consistent with the experimental results obtained by means of PIV. Figure 13 - Figure 16 present the time-averaged flow fields obtained with the use of CFD for the analyzed methods and mesh types zoomed in the region of the contact point between the first and second granule. Velocity magnitude as well as stream-wise and span-wise velocities are depicted in order to compare different cases of CFD and PIV results in detail. It can be noticed that the influence of granules' contact point representation method can be seen only in the region occupied by the additional cylindrical volume generated with the CPE method where the velocity does not exceed 0.02m/s (less than 1% of inflow velocity). However, variations can be observed regarding the mesh type: higher velocity values (more similar to PIV) are observed in the case of TETRA, but POLY represents the boundary layer with better agreement of the experimental research results. In general, CFD underrates velocity in the region between the granules in comparison to PIV. The most significant differences among different cases of CFD simulations are observed between GRR POLY and GRR TETRA, which proves the importance of the applied mesh type.

Conclusions
The presented results prove that the implementation of the CPE method does not affect the flow above the granules. A reduction of granule diameter in the GRR method slightly affects the flow above the granules but has significant effect on the calculated porosity of the packed bed of granular material which increases 3.26%. In contrast, the CPE method affects the calculated porosity only by 0.03%. The local velocity in the region close to the granules' contact point is underrated by CFD in comparison to PIV. The applied mesh type influences the flow field more significantly than the applied method of contact point representation. The better characteristics of POLY mesh in terms of numerical diffusivity is revealed in more accurate boundary layer representation and better agreement of CFD and PIV results regarding velocity fluctuations particularly along the vertical axis of the domain. Moreover the results are generated within approx. 2.5 times lower number of iterations in the case of POLY mesh, which may be a major advantage for industrial applications.