Computational domain discretization in numerical analysis of flow within granular materials

The discretization of computational domain is a crucial step in Computational Fluid Dynamics (CFD) because it influences not only the numerical stability of the analysed model but also the agreement of obtained results and real data. Modelling flow in packed beds of granular materials is a very challenging task in terms of discretization due to the existence of narrow spaces between spherical granules contacting tangentially in a single point. Standard approach to this issue results in a low quality mesh and unreliable results in consequence. Therefore the common method is to reduce the diameter of the modelled granules in order to eliminate the single-point contact between the individual granules. The drawback of such method is the adulteration of flow and contact heat resistance among others. Therefore an innovative method is proposed in the paper: single-point contact is extended to a cylinder-shaped volume contact. Such approach eliminates the low quality mesh elements and simultaneously introduces only slight distortion to the flow as well as contact heat transfer. The performed analysis of numerous test cases prove the great potential of the proposed method of meshing the packed beds of granular materials.


Introduction
The widespread prevalence of granular materials in nature and their common use in many industries, make packed beds of granular materials the focus of numerous experimental and theoretical studies [1][2][3][4]. Such research are often performed using advanced computational techniques [5][6][7], including Computational Fluid Dynamics (CFD), which applies numerical methods to solve fluid flow problems including Fluid-Solid Interactions (FSI).
Numerical methods of solving the differential equations are implemented in commercial CFD codes but the user has been left with detailed control over the discretization of the computational domain. This process is commonly called meshing and consists in dividing the analysed geometry into numerous small control volumes commonly called cells or elements. This step of CFD influences the numerical stability of the model and consequently the accuracy of the obtained results [8,9].
Moreover the common challenge for CFD is the expensive computational cost in terms of time required to reach a converged solution [10]. The very first issue influencing the convergence as well as the accuracy of the results is the applied method of computational domain discretization [8,11,12].
Therefore the carried out research focus on the innovative approach to discretization method of the analysed domain representing the packed bed of granular material.

Methods
The performed research concern the development of the optimal method of discretization of the computational domain which represents solid granules submerged in the fluid domain -packed bed of granular material. The solid bodies of spherical shape (homogeneous granules) are in direct contact between each other in a single point between two adjacent granules.
The standard approach to computational domain discretization where two spherical bodies contact tangentially in one point (Fig. 1a) leads to extremely low-quality mesh cells, resulting in numerical instability of the flow calculations. The reduction of radius of individual granules (Fig.  1b) can be applied in order to avoid the above mentioned effect of low quality mesh. It allows to obtain the acceptable mesh quality but entails the necessity of significant size reduction of mesh cells. In consequence the total number of cells increases and it directly affects the overall calculation time. Unfortunately this approach disturbs the flow and heat transfer between individual granules which no longer have direct contact with each other.
Therefore an innovative approach to discretization of packed beds of granular materials is proposed. It consists in extending the contact between the individual granules from single point to the cylindrical region as shown in Fig.  1c. Such approach allows to achieve high quality mesh and simultaneously guarantees the contact of the adjacent granules. Therefore the flow disturbance between the granules is insignificant.

Geometry and analysed cases
The geometry representing the packed bed of granular material ( Fig. 2) was generated using SolidWorks with ANSYS 18.1 plugin, which allows the bi-directional associativity between the CAD geometry (SolidWorks) and mesh generator (ANSYS Meshing). Two global variables defining the relative radius of individual granule (eq. 1) as well as relative radius of the cylindrical volume (eq. 2) were introduced in order to speed up the process of preparing the geometry for different cases.
where : R radius of the granules contacting in a single point; Rg reduced radius of individual granules; Rcyl radius of the additional cylindrical volume between individual granules.
The geometry was prepared in two configurations in order to generate the mesh using two methods ( Fig. 3): -with additional cylindrical volumes between the individual granules, -with reduction of granules radius. Five cases (a-e) were analysed for each method according to Table 1. The CAD model of the packed bed was prepared as a sector of the whole bed consisting of one spherical granule (marked red in Fig. 3) surrounded by eight sections of granules. Each section represents one-eight of a granule (marked green in Fig. 3). Such procedure allows to minimize the volume of the computational domain which directly results in CPU time reduction. Simultaneously the possibility to generate the entire packed bed by simple multiplication of the sector CAD model is maintained.

Meshing
The meshing process of all ten configurations of the analysed geometries (five cases for each of the two methods) was performed using ANSYS Meshing 18.1. Tetrahedral (TET) elements were used to discretize the computational domain. The size function was set to UNIFORM in order to obtain the most regular mesh. Moreover the minimal size, maximal face size as well as maximal TET size constrains were set to the value of 0.1R. The above meshing settings resulted in the final mesh consisting of 36k to 47k elements (depending on the method and case). The final mesh for selected case of each of the two methods is depicted in Fig. 4.
The output mesh quality was analysed. The quality analysis provides a composite metric that ranges between 0 and 1 according to formula (3). A value of 1 indicates a perfect element while a value of 0 indicates that the element has a zero or negative volume.
The element of the worst quality generated throughout the entire research (ten meshes) was characterized by the quality index of 0.17 and therefore it was assumed that the mesh did not affect the numerical stability of the computational model.
where : V volume of the computational domain; l length of the element edge; C constant dependent on the element type (eg. 124.708 for TETs).

Computational model
The exactly identical computational model (except for the domain geometry) was used for all analysed cases and both discretization methods. The solver was configured as pressure-based and the analysis were performed for steady state. 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.
Air at normal conditions was the fluid medium. Velocity-inlet boundary condition type was assigned to the inlet (magenta surface in Fig. 5) with velocity magnitude normal to the boundary equal 1 m/s. The outlet of the computational domain (yellow surface in Fig. 5) was defined as outflow. Wall with no slip shear condition was assigned to the surfaces of granules (red and green surfaces in Fig. 5) whereas wall with specified shear condition of shear stress X, Y and Z components equal 0 was assigned to the surfaces of module interfaces (cyan surfaces in Fig. 5) in order to assure the periodicity of the flow in each of the two directions. The calculations were performed for only 100 iterations in order to compare the rate of convergence for all analysed cases.

Results
The obtained results concerning the velocity magnitude averaged along outlet line and averaged on outlet surface are depicted in Fig. 6 and Fig. 7 respectively. Cases a-e for both methods are parametrized according to Table 1. The average velocity magnitude is referenced to the inlet velocity.  The discretization method with additional cylindrical volume is characterized by relatively stable average velocity along outlet line (variation of 3%) and on outlet surface (variation of 1%) for each of the analysed radius of the cylindrical volume (Rcyl%). The variations of average velocity are more significant (18% and 5% respectively) for the method of granule radius reduction.
Contour plots of velocity magnitude on outlet surface, XY and XZ planes for both methods and all configurations are presented in Fig. 8-Fig. 13.  Fig. 8. Velocity magnitude on outlet surface obtained with the discretization method of granule radius reduction Rg% equal: a -99.75%. b -99.50%. c -99.25%. d -99.00%. e -95.00%.   It is clearly seen that the results obtained with the mesh prepared according to the proposed innovative method of computational domain discretization are similar irrespective of the radius of the additional cylindrical volume (Rcyl%) for all post-processing surfaces.
The method of granule radius reduction is sensitive to the level of reduction Rg% -the velocity field becomes more uniform along with the decrease of granule radius.

Discussion
The proposed innovative method of computational domain discretisation in CFD analysis of packed beds of granular materials is characterized by lower variation of mesh elements count (37 325 -36 590) in comparison to the commonly used method of granules radius reduction (37 731 -47 300).
The variations of average velocity for the analysed cases are significant for the granules radius reduction method and negligible for the innovative method proposed in the paper.
The same issue occurs in terms of the velocity field obtained using the mesh generated with the analysed methods: the proposed method is independent to the radius of the additional cylindrical volume (Rcyl%) while the level of reduction of granule radius (Rg%) in the standard method strongly affects the velocity field.
Experimental validation of the numerical results is necessary in order to properly and fully define the applicability of the presented method. The experimental test stand has already been designed and it is under construction.

Conclusions
It can be concluded that the proposed innovative method of computational domain discretization is characterized by more stable behaviour of the CFD results in comparison to the commonly used method.
The carried out research indicate the purposefulness of further investigation in order to evaluate the optimal Rcyl% value as well as the correlation of mesh element size and the radius of the additional cylindrical volume.
Moreover the Grid Convergence Index (GCI) [13] has to be checked to assure the mesh-independent solution and fully evaluate the proposed method.
Finally the experimental research is necessary to validate the method against real data.