Simulation of Heat and Mass Transfer in a Grain Pile on the Basis of a 2D Irregular Pore Network

: The so-called pore network model has great advantages in describing the process of heat and mass transfer in porous media. In order to construct a random two-dimensional (2D) irregular pore network model for an unconsolidated material, image processing technology was used to extract the required topological and geometric information from a 2D sample of soybean particles, and a dedicated algorithm was elaborated to merge some adjacent small pores. Based on the extracted information, a 2D pore network model including particle information was reconstructed and verified to reflect the pore structure of discrete particles. This method was used to reconstruct a random 2D irregular pore network model of wheat. Accordingly, a multi-scale heat and mass transfer model was implemented to simulate the drying of wheat. The simulation results were consistent with the experimental results, which indicates that the reconstructed irregular pore network model can effectively simulate the real pore structure inside unconsolidated porous media. The present approach may be regarded as the foundation for establishing in the future a three-dimensional pore network model and studying the heat and mass transfer process in a grain pile.


Nomenclature
Wt S S1 SP Lt y l ε width of the throat rectangular area total area of particles in 2D image total area of pores total length of throats the number of throats at a certain length length of a throat, m porosity wet air density, kg/m 3 permeability, m 2 wet air viscosity, Pa•s laplacian operator wet air pressure, Pa throat radius, m total number of pores connected to the pore i1 volume of a pore, m 3 vapor density within pore, kg/m 3 air speed in throat, m/s surface area of particle, m 2 phenomenological coefficient particle radius, m the thickness of the single layer after dividing the particle into mt layers along the radius jth calculation on the particle scale the interval between each calculation time, ∆t=10Δτ convective mass transfer coefficient, m/s vapor density at particle surface, kg/m 3 total number of particle around the pore i3 number of pores around a particle sign of wind flowing into a pore sign of wind out of a pore air temperature in a pore, K specific heat capacity of air, kJ/(kg•K) air density, kg/m 3 convective heat transfer coefficient, W/(m 2 •K) particle temperature, K mass diffusion coefficient, m 2 /s spherical coordinates dry-basis moisture content of particle dry-basis equilibrium moisture content of particle latent heat of water, kJ/kg specific heat capacity of vapor, kJ/(kg•K) mass of evaporated water from grain particle in unit time, kg/s specific heat capacity of water, kJ/(kg•K) specific heat capacity of dry grain particle, kJ/(kg•K) Cg ρg density of dry grain particle, kg/m 3 of porous media. In the early 1990s, Daian and Saliba introduced the pore network method into the research field of porous media drying [Daian and Saliba (1991)]. In the study of the drying process of rigid porous media, Prat [Prat (1993)] combined seepage theory with pore network theory and obtained the phase distribution and drying rate curve at different times during the drying of porous media. The combination of pore network theory and drying theory has promoted the application of pore network theory in simulating the drying process [Laurindo and Prat (1998); Laurindo and Prat (1996); Prat (1995); Segura and Toledo (2005)]. In 2006, Yang [Yang (2006)] established a mathematical drying model based on a 2D fractal pore network, and the simulated drying curve of a potato slice was more consistent with the experimental results than the drying model based on the regular pore network. Gong et al. [Gong, Niu, Dong et al. (2009)] obtained the microstructure parameters of the frozen-dried potato through a mercury intrusion experiment, and the fractal dimension of pores by calculation. The image reconstruction pore network model is completely based on the scanned threedimensional images or the combination of a series two-dimensional images. The specific locations of pores and throats should be determined during reconstruction. Therefore, it is the closest to the real pore structure. It is suitable for developing a small pore network model limited by the image size. Based on actual solid particle accumulation, Wang et al. [Wang, Kharaghani, Metzger et al. (2012)] obtained a three-dimensional X-ray scanned image, extracted the pore network information, and simulated the liquid transfer and drying processes in pores. The simulation results were largely consistent with the experimental results. However, the pore network information comes from the actual stacking material and cannot be used to reconstruct a pore network model with a random size. To solve this problem, Zhao et al. [Zhao, Liu, Han et al. (2018)] developed a threedimensional reconstruction technique based on the marching cubes algorithm and ray casting algorithm, which restored the geometry and spatial distribution of soil pores. The ray casting algorithm can clearly reflect the pore contour and the detailed pore structure information compared with the marching cubes algorithm. Most of the processes of heat and mass transfer only occurred in pore space but not between the skeleton and pore space in the above studies, which is not in line with the actual situation, especially in the hygroscopic porous media. Therefore, it is necessary to combine the particle information with the pore and throat information in the porous medium. Yuan et al. [Yuan, Tan, Zhao et al. (2016)] constructed a regular pore network of unconsolidated corn kernels by statistical methods and established the multiscale heat and mass transfer models for simulating the drying process of corn. Yue et al. [Yue and Zhang (2017)] developed an algorithm to construct a tetrahedral model, equivalently obtained the pore space and pore structures in bulk grain, and then calculated the airflow distribution in the pore space. Kharaghani et al. ] constructed an irregular pore network model for convective drying of random aggregates composed of mono-sized primary particles, and an isothermal pore network drying model based on invasion percolation driven by evaporation. The simulation results are in line with the actual drying conditions. However, the moisture is dried only in pore space and does not exchange with the skeleton. There are some shorter throats in the pore network model, which are not merged and inconsistent with the actual structure.
In this study, soybean was selected as the test material for extracting pore space information and particle information from a 2D unconsolidated porous medium by MATLAB image processing technology. Based on Thiessen polygons, an irregular 2D pore network model reflecting the real accumulation structure of the grain pile was constructed, which can be seen as a pore network structure used to simulate any nonconsolidated porous media stacked structure. Combined the constructed irregular pore network model, a heat and mass transfer model was established to simulate the drying process of particulate material and verified by a wheat drying experiment. It is helpful for constructing a 2D and 3D irregular pore network of unconsolidated material, simulating the heat and mass transfer process in porous media.
2 Construction of a 2D irregular pore network

Extraction of structural parameters from an 2D image
Most agricultural materials are irregular and the size of an irregularly shaped object can be described by the tri-axial dimensions-length (L), width (W) and thickness (T). It is difficult to address irregular particles, so an irregular particle is often taken as a sphere, and the equivalent radius Rp can be calculated by Eq. (1).
Monolayer soybeans were randomly placed in a 0.1 m×0.1 m area until soybeans filled the entire rectangular space. The scanned 2D image of monolayer soybeans was processed by MATLAB software. Through binary image processing Fig. 1(b), particle skeleton segmentation Fig. 1(c), and watershed segmentation Fig. 1(d), the particle information and the topological structure of the pore space in the 2D image were obtained. The particle number, position and area can be determined by particle skeleton segmentation and the connected region tag function, and the center and equivalent diameter of an equivalent circle with the same area of the particle can be obtained. The intersections in Fig. 1(d) are considered as the centers of pores, and the edges are considered as the centerline between two adjacent pores. During watershed segmentation, some short edges can appear, as shown in Fig. 2(a). The vertices A1 and A2 represent the centers of two adjacent pores, but the two small pores should belong to the same larger pore. Therefore, it is necessary to identify this situation and merge the small pores to obtain the pore parameters that are closer to the actual pore structure. The midpoint A3 of A1, A2 is the new center of the merged pore, as shown in Fig. 2 The minimum pore appears when the three circles surrounding the pore are tangent to each other, as shown in Fig. 3, and the minimum coordination number is 3. The minimum pore radius r can be obtained by Eq.
(3) Figure 3: Schematic diagram of a minimum pore Therefore, the edges smaller than 2r cannot reflect the actual throat condition, and the small pores should be merged and a new center can be obtained. The pore radius can be extracted based on the maximum sphere algorithm, and the throat length can be calculated though the centerline spacing and the radii of two adjacent pores.
Assuming the width of all throats is consistent, and the throat width can be calculated by Eq. (4). The coordination number is equal to the number of edges connected to a pore.
According to the above processing, 193 particles, 382 pores, and 573 throats were obtained from Fig. 1. The average coordination number is 3.00. The equivalent radius distributions of soybean particles and extracted particles from the image are shown in Fig.  4. The average equivalent radius of extracted particles is 3.4443 mm, which is consistent with the average equivalent radius of soybean particles of 3.4496 mm.

Figure 4: Equivalent radius distribution of soybeans
The extracted parameter distribution of the coordination number, the throat length and the pore radius are shown in Figs. 7-9 respectively. The throat length mainly distributed between 3 mm and 7 mm with the normal distribution. After fitting, the throat length distribution roughly conforms to formula (5) The pore radius generally varies between 0.7 mm and 1.5 mm, as shown in Fig. 8, which is smaller than the radius of soybeans. There are also some large pores caused by friction, as the soybean surface is not completely smooth and the shape is not completely spherical. The calculated porosity from the image is 27.85%, and the throat width is 0.2761 mm according to Eq. (4). The extracted porosity is smaller than the real porosity of the grain pile because the pore space in the 2D image is just the minimum projection area of the real pore space on the normal projection plane. Therefore, the 2D pore network model is limited in simulating the three-dimensional pore space.

Construction of a 2D irregular pore network model based on the extracted parameters
By analyzing the extracted parameters and considering the narrow distribution of the equivalent diameter and the effect of image noise on the extracted parameters, it follows that soybean particles can be regarded as equal-diameter spheres with a radius of 3.44 mm. According to the principle of Thiessen polygons, this paper argues that this construction can accurately describe the 2D soybean accumulation. OpenGL technology is used to sprinkle a certain number of equal-diameter circles randomly inside the fixed rectangular area. Moreover, all of the equal-diameter circles do not intersect each other. Therefore, the position of each circle center can be obtained, and the central distance between two adjacent circles is greater than or equal to the diameter 2R of the circle. Based on above circle centers, Thiessen polygons are constructed, and the topological information and geometric information are obtained following the similar process of extraction. A new value is identified and reassigned to the short edges with the length less than 2r according to Eq. (3), and the topology information and geometric information of the 2D irregular pore network model are updated. The flow chart of the construction program is shown in Fig. 5.  6 shows the constructed 2D irregular pore network model in a 0.1 m×0.1 m area. The vertex and edge of a polygon represent the pore and throat, respectively, and the closed space in the polygon represents the particle region. Comparing the extracted parameters from a 2D image and the parameters in the reconstructed 2D irregular pore network, all the parameter distributions are very similar, as shown in Figs. 6-8. The reconstructed coordination number and the throat length are consistent with the extracted parameters, and some large pores appear in the reconstructed pore network model caused by the uneven distribution near the rectangular edge during the reconstruction process. The differences between the extracted parameters and the reconstructed parameters are shown in Tab. 1.     There are 193 particle regions, 388 pores and 580 throats in the reconstructed 2D irregular pore network model. In addition, the porosity of the pore network model is 28.25%. The throat width is 0.2661 mm, which is slightly larger than the extracted parameters. The average coordination number is 2.99, which is slightly smaller than the extracted parameters. The maximum error of the throat width is 3.62%, and the minimum error of the average coordination number is 0.33%, while the errors of other parameters are less than 2%. Therefore, the reconstructed model can well reflect the pore structure of granular materials and can be used to construct a large-scale pore network model. The previous 2D pore network models based on a slice were mostly regular pore network models. The particle radius, pore radius and throat length are assigned according to their distributions and a certain rule, which is considered equivalent to an irregular pore network model. However, the particle radius distribution obtained from a slice is very different from the real particle radius distribution, and the coordination number of a regular pore network model is fixed, which is 4 in a frequently used regular rectangular pore network model except for the boundary. It cannot match the particle size and the actual pore network parameters in bulk grain. The method of constructing a 2D pore network model in this paper is different from the previous method, and it can reflect the real particle size and the pore network parameters in a 2D image.

Construction of a 2D irregular pore network model of wheat
A wheat kernel is approximately ellipsoid; its equivalent radius of more than 93.5% of wheat particles is relatively concentrated at 1.90~2.20 mm in this study, as shown in Fig.  10. Although there is a certain deviation between the shape of a wheat kernel and a sphere, wheat particles are taken as equal-diameter spheres with an equivalent radius of 2 mm to construct the pore network model. Based on the above method and the equivalent size of a wheat kernel, a 2D irregular pore network model of wheat is constructed in a 0.1 m×0.1 m container, as shown in Fig. 11. Therefore, 631 wheat kernels, 1264 pores and 1894 throats are obtained, and the average coordination number is 2.997.

Simulation and experimental study of wheat drying
Based on the 2D irregular pore network model of wheat, the heat and mass transfer process during wheat drying is explored by simulation and experiments. A multiscale mathematical model is established to describe the wheat drying process, which is shown in Eqs. (7), (8), (10) and (11). A drying experiment is carried out to verify the simulation results.

Multiscale drying model based on the pore network model
Assumptions: (1) Container scale a) Heat transfer is dominated by convection, conduction and radiation are ignored. b) Mass transfer is dominated by convection, diffusion and micro-scale effects are ignored. c) Airflow and vapor are represented as incompressible ideal gas. d) Heat and mass transfer are insulated by the container wall.
(2) Particle scale a) No grain particle shrinkage and deformation during drying. b) Uniform temperature distribution in grain particle during drying. c) Moisture distribution in a grain particle depends on liquid or vapor diffusion.

Momentum transfer equation
Darcy's law is used to describe the airflow passing through grain pile.
The momentum equation can be obtained by discretizing the above equation:

Heat and mass transfer equation on particle scale
Since the evaporation rate from the grain surface is greater than the moisture diffusion rate in the grain, the grain drying rate decreases, and the grain temperature increases until it reaches the air temperature surrounding the grain. Luikov proposed the following equations to describe the slow drying process [Cao and Zhu (2001)].
The value of G is related to temperature, moisture and pressure. To simplify the calculation, the grain particle is taken as continuous media, as shown in Fig. 13. Combined the heat conservation formula, the discrete heat and mass transfer equation on the particle scale is as follows.   Fig. 14 shows the heat and mass exchange between the two scales; and the heat and mass exchange control equation can be obtained.

Heat and mass exchange equation between two scales
These equations can be simplified into series of linear equations and solved by MATLAB software.

Wheat drying experiment 3.2.1 Materials, instruments and equipment
Material: The initial moisture content of wheat is too low to be dried. Therefore, the initial moisture content should be raised by humidification. Wheat is immersed into water for 24 hours, spread on the floor into a 2-cm thickness and turned frequently until its drybasis moisture content is approximately 25%, and then piled to make it uniform. Instruments: Electronic balance (FB224, Shunyuhengping, China), hand-held temperature and humidity meter (HMP76, Vaisala, Finland), and infrared thermal imager (Testo 875-1i, Testo, Germany). Equipment: Oven (DHG-9140A, Jinghong, China), and small bin with the length, width and height of 0.6 m×0.6 m×1.2 m as shown in Fig. 15.

Operating conditions and methods
Drying experiments were carried out in the bin (Fig. 15) under the operation conditions shown in Tab. 2. A U-shape orifice plate was used to simulate the U-shape ventilation duct structure in real bins. The grain temperature and moisture content at different sampling positions and time intervals were measured quickly to minimize the effect of environmental conditions. The sampling device and 4 sampling positions at the same grain layer are shown in Fig. 16. The holes in the sampler determine the sampling positions in the depth direction.

Experimental results and analysis
Limited by the computational workload, only the material with a layer thickness of 10 cm close to the orifice is selected as the simulation area in this study. Fig. 17 shows the moisture content of wheat at different drying times and sampling positions at layer D.
The moisture content of wheat at the same layer is almost same, and the unevenness of the moisture content is less than 2.5%. Comparing the drying curves at 4 sampling points, the drying rate at point 2 is larger than that at point 1. It means that the drying rate at the inlet boundary is less than that at the inlet corner. The wheat temperatures at layer D are almost consistent and change with drying time, as shown in Fig. 18. The wheat temperature rises rapidly at first, and then slowly becomes close to the ambient air temperature. The wheat temperature at point 2 is obviously higher than that at point 4.

Simulation of wheat drying process 3.3.1 Initial conditions and boundary conditions
The initial dry-basis moisture content of wheat was 0.2353; the initial wheat temperature was 293.81 K; the inlet air pressure was 52.18 Pa, the vapor density was 0.00169 kg/m 3 ; the inlet air temperature was 335.15 K; and the initial internal air temperature was 293.81 K. The center coordinates of the sampled pores and particles are listed in Tab. 3.  Fig. 19 compares the simulated average moisture content with the experimental results. The simulated moisture content decreased from 0.2353 to 0.0710 in approximately 9.3 hours, which is slightly faster than the experiment. In the early drying stage, the simulated drying rate is similar to the experimental drying rate. However, the simulated drying rate slows down as it approaches a later stage. When the final moisture content is the same, the drying time required for the experiment and simulation is almost the same.

Analysis of the simulation results
The experimental drying rate that is faster than the simulated results may be caused by the sampling process, which will disturb the stable state and form an airflow channel with a larger porosity than other areas in the grain pile. Fig. 20 depicts the vapor densities in different sampling points in the three-layer depth. The closer to the bottom, the larger the vapor density is. This is because the vapor at the top layer is easily removed from the thin grain pile by hot air. All vapor densities in different sampling points initially increase with the moisture being removed easily and quickly from the wheat at the beginning stage, but gradually decrease when the moisture content decreases and the water inside the wheat becomes increasingly more difficult to be removed at the later drying stage. The closer to the bottom, the higher the wheat temperature is, which occurs because wheat at the low layer is heated by the hottest air. The average simulated wheat temperature is higher than the average wheat temperature measured at 4 sampling points, as shown in Fig. 16. The measured low wheat temperature might be caused by the instrument test accuracy (±2 o C), measurement error (wheat temperature decreases quickly under room temperature in winter), and the initial temperature error between the actual wheat temperature and the assigned value. The wheat temperature gradually increases to the inlet air temperature. The air temperatures in the sampling pores can be observed in Fig. 22. As a large amount of heat required for vaporization is absorbed from air, the air temperature decreases at the beginning drying stage, and then increases and approaches the inlet air temperature with less moisture evaporation.

Figure 22: Air temperatures in sampling pores
4 Conclusions 1) A 2D irregular pore network model was constructed based on the image of monolayer material, and the distribution of extracted particle sizes is similar to that of the real particles.
2) The pore space parameters of the reconstructed irregular pore network model are similar to the extracted parameters from a 2D image. The reconstruction method can be used for reconstruction of an irregular pore network model with a random size.
3) The established multiscale heat and mass transfer model can simulate the drying process of wheat and be verified by experiments. It proves that the pore network model can effectively describe the actual pore structure of stacking unconsolidated material.