Changes in Pore Structures of Porous Beds When Subjected to Vertical Vibration †

The objective of this study was to investigate the effect of vibration on critical pore structure parameters related to flow in porous beds. The discrete element method was used to simulate particle packing in porous beds of soybean subjected to vibration. The porous bed was simulated as an assembly of spherical particles with diameters randomly distributed between 5.5 mm and 7.5 mm. The simulated porous bed was subjected to vertical vibration at a fixed frequency of 15 Hz and multiple amplitudes from 0.5 to 4.0 mm, resulting in vibration intensities from 0.45g to 3.62g (g = gravitational acceleration). The location (coordinates) of each particle was tracked during vibration. Based on the simulated spatial arrangement of particles, critical flow-related parameters of the porous bed, including porosity, tortuosity, and pore throat width were calculated. It was found that vibration intensity of 1.81g resulted in the lowest porosity, whereas lower vibration intensity did not have enough energy to densify the bed and higher intensity produced less dense parking due to over-excitement. Local porosity fluctuated markedly during vibration, with a general trend of decrease as vibration progressed. Vibration noticeably affected the shape (tortuousness) of flow path. Tortuosity of the porous bed before vibration was higher (2 % to 9 %) than that after vibration. Vibration reduced the pore throat width by 18 % on average (from 3.3 mm before vibration to an average of 2.7 mm after vibration).


Introduction
The knowledge of pore structure of porous beds (media) is critical in understanding fluid flow through porous beds.A number of 2D (Lao H.W. et al., 2004;Ljung A.L. et al., 2012;Xiao Z.F. et al., 2008;Yiotis A.G. et al., 2010) and 3D models (Balhoff M.T. et al., 2007;Sobieski W. et al., 2012;Thompson K.E., 2002;Thompson K.E. and Fogler H.S., 1997) have been developed to describe pore structures of porous media.Those models generally deal with porous beds with fixed pore structures.However, porous beds are often subject to such mechanical disturbances as vibration during operation.When a porous bed consisting of randomly packed particles is vibrated, the pore structure is altered due to re-arrangement of particles.The mechanisms by which the pore structure is affected by vibration are extremely complex.An X. Z. et al. (2005) used the discrete element method to study the densification of porous beds (sphere packing) when subjected to vertical vibration.They showed that the degree of densification was significantly affected by both vibration frequency and amplitude, and identified two densification mechanisms: pushing filling and jumping filling.Pushing filling generally occurred at low vibration intensity in which the contact between spheres was maintained, while jumping filling at high intensity in which the contact between particles was periodically broken.An X. Z. et al. (2009) conducted experiments to study the effect of vibration on packing density and observed that there existed optimum values of vibration frequency and amplitude to achieve the maximum packing density.Frequency and amplitude should be considered separately in terms of their effect on packing density.Liu C. et al. (2010) used a high speed camera to study the behavior of particles in a bin subjected to vertical vibration.They reported that all particles in a porous bed vibrated as a whole (a continuum) when the bed was subjected to low intensity vibration, whereas individual particles vibrated in their own modes (as a discrete system) when subjected to high intensity vibration.
The effect of vibration on pore structures is complex.Most studies on vibration of porous beds have been focused on particle packing, and few have explored how changes in particle packing caused by vibration would impact such properties as porosity and tortuosity which are related to flow through the porous bed.Traditionally, the studies of particle packing are pursued in the realm of particulate mechanics, while flow through porous media has established itself as a sub-discipline in fluid mechanics.This study attempted to establish a connection between flow in porous media and particulate mechanics.A powerful tool used in particulate mechanics-the discrete element method, was employed to study pore structures from the angle of flow in porous media.The objective of this paper was to investigate the effect of vibration on pore structures of porous beds with an emphasis on parameters that are related to flow through the porous beds.Specifically, a discrete element model (DEM) was used to simulate the formation of porous beds and changes in pore structure during vibration; and an algorithm was developed to use the DEM simulation data to quantify the flow-related characteristics of the porous bed, including porosity, flow path (length and shape), tortuosity, and pore throat diameter.

Discrete element simulations of porous bed subjected to vibration
The discrete element method was originally proposed by Cundall P.A. and Strack O.D.L. (1979) to solve problems in rock mechanics.The method has been widely used to analyze granular and discontinuous materials.A commercial discrete element software package PFC 3D (Version 4.0 EV, Itasca Consulting Group Inc., Minneapolis, USA) was used to build a simulation model in this study.The PFC 3D model first simulated the formation of a porous bed, during which the trajectory and rotation of each particle were tracked to evaluate its position and orientation in a time stepping simulation.Once the porous bed was formed, a vibration excitation was applied.Simulations were conducted for a model bin 0.28 m in height and 0.15 m in diameter, filled with soybean to a depth of 0.25 m (Fig. 1).This model bin was used by Liu C. et al. (2010) to study the vibration of soybean in storage.The same mechanical property parameters for the PFC 3D model (including the stored soybean) used by Liu C. et al. (2010) were adopted in this study (Table 1).Based on the measured sizes of bin and soybean kernels, the model bin was simulated as a cylindrical container filled with spherical particles.It should be noted that soybean kernels are not spherical, but they resemble a sphere with high average sphericity values (above 0.8), and therefore, most researchers have used spheres to approximate soybean kernels when simulating bulk soybean characteristics in DEM models to reduce computation times (Boac J.M. et al., 2014).Based on the measurements of 24 randomly selected kernels, the average long, intermediate, and short axis dimensions for soybean used in this study were 7.6, 6.7, and 5.7 mm, respectively.The corresponding geometric mean diameter was 6.6 mm, and sphericity 0.87.Considering variations in measured soybean kernel size, a rage of particle sizes from 5.5 and 7.5 mm recommended by Liu C. et al. (2010) for DEM simulations was adopted in this study.A coordinate system was set up with x and y representing the horizontal directions, z the vertical direction, and the origin at the bottom-center (Fig. 1).
A simple harmonic excitation was used to vibrate the bin.Vibration excitation was applied to the bin floor (bottom) in the vertical direction with prescribed frequencies and amplitudes as follows: (1)  A frequency of 15 Hz and amplitudes of 0.5, 1.0, 2.0, 3.0, and 4.0 mm were used in simulations.These frequency and amplitude combinations were selected to achieve the maximum densification.According to the experimental results reported by An X. Z. et al. (2009), the dense packing was achieved with vibration intensity between 1.5g and 4.0g (the peak density occurred at 2.0g-2.5g).The vibration intensity levels, determined as A(2πf ) 2 , were 0.45g, 0.91g, 1.81g, 2.72g, and 3.62g for amplitudes of 0.5, 1.0, 2.0, 3.0, and 4.0 mm at 15 Hz, respectively.
Vibration was simulated in time steps of 3.3 μs.This time step (Δt) was determined from the mass (m) and the stiffness (k pn ) of particles (Δt = pn / m k ) (Itasca, 2005).Preliminary simulations were conducted to determine the total length of vibration time required to reach the maximum density and it was found the 16 s was sufficient for all vibration intensities.At each time step, the position (x, y and z coordinates) of each particle was saved in a data file for late use in calculating flow parameters.
One of the most important parameters affecting flow through a porous bed is porosity.In this study, both global (overall) and local porosity values were calculated based on DEM simulation data.The global porosity was calculated as the ratio of the total volume occupied by the particles to the total volume of porous bed as follows: In a porous bed consisting of randomly packed particles, the pore distribution would not be uniform.Therefore, local porosity was also determined using the built-in function in PFC 3D at five locations, which were denoted as B1(0.06,0, 0.119), B2(0.06, 0, 0.129), B3(0.06, 0, 0.130), B4(0.06, 0, 0.136), and B5(0.06, 0, 0.146) (three numbers in parenthesis are x-, y-and z-coordinates, respectively) (Fig. 1).Because the bin was relatively small, these five locations were selected near the bin center to avoid the boundary effects.The instantaneous local porosity values at these locations were calculated at every time step (3.3 μs) during vibration simulation.

Airflow path
The airflow path was calculated using the methodology proposed by Yue R. and Zhang Q. (2013).The porous bed was treated as a collection of tetrahedron units, and each tetrahedron contained a pore surrounded by four particles (Fig. 2a).Air (fluid) was assumed to enter the tetrahedron unit from the base (inlet) triangle (ΔABC) and exit from the other three (exit) triangles.The lines connecting the centroid of the base triangle to the centroids of the three exit triangles represented three local airflow paths.The length of each local path was calculated and the shortest and longest paths were identified.To compare the width of each airflow path, the open area (pore throat) between the three particles that form a triangular face in the tetrahedron unit was calculated for all four face triangles (Fig. 2b).The open area was calculated as the area of face triangle less the projected area of three particles on the triangle face (Dudda W. and Sobieski W., 2014).The equivalent pore throat width was considered as the diameter of a circle which had the same area as the open projected area.
Based on the calculated pore throat width, the widest and the narrowest local flow paths were identified.To construct a specific type of global flow path (running from the bottom to the top of the porous bed), say the widest path, one of the three exit triangles that had the widest pore throat was selected as the base (inlet) triangle for the next tetrahedron unit.This process was repeated from the bottom to the top surface of the porous bed and a global flow path was obtained by connecting the local paths in all associated tetrahedron units (Fig. 3).It should be noted that simply connecting a series of straight lines (local paths) to form a global path resulted in sharp turn angles between local paths (Fig. 3).However, flowing air cannot make sharp turns.Therefore, a procedure proposed by Sobieski W. et al. (2012) was used to replace each sharp turn angle with an arc (Fig. 3).

Tortuosity
Different definitions of tortuosity can be found in the literature (Bear J., 2013).The geometric tortuosity was used in this study, which was defined as the ratio of the actual flow path length to the depth of porous bed: The actual length of a (global) flow path was determined as the sum of all local path lengths.The depth of porous bed was determined from the positions of particles in the top layer of bed simulated by the DEM model as follows (Sobieski W. et al., 2012):

Experiment
A cylindrical model bin of 0.28 m height and 0.15 m diameter was used to conduct tests to verify the DEM simulations (Fig. 4).The bin was made of Plexiglas and placed on a vibration table (Model 307-97056, Cougar Industries Inc., Peru, IL, USA).The table was supported by three piston-spring assemblies to permit only vertical displacement.The vibration table was driven by a pneumatic piston with compressed air and controlled by a pressure regulator to achieve different vibration intensities.A portable vibration analyzer (Model 654 CS, Vitec Inc., Cleveland, OH, USA) was used to measure the vibration amplitude and frequency and the signal was displayed on an oscilloscope (Model 190 B/C ScopeMeter, Fluke Industrial B. V., the Netherlands).
For each test, measured amount (mass) of soybean was placed into the bin through a funnel located 350 mm above the bin, and the top surface was levelled manually after filling.The total volume of the porous bed was calculated from the filling depth and the bin diameter, and the packing density (ρ b ) was then calculated from the total volume and the mass of soybean in the bin.All linear dimensions were measured by using a ruler with a resolution of 1 mm, and mass measured by an electronic scale with a resolution of 0.1g.
The liquid displacement method was used to measure the particle density (ρ p ) of soybean kernels.The specific steps were as follows: (1) 100 kernels were weighed and placed into a 50-mL flask filled with toluene, (2) the volume of 100 kernels was determined from the displaced volume of toluene, and (3) the particle density was calculated from the measured kernel mass and volume.A steel calibration ball of 8580 mm 3 (25.4mm in diameter) was used to calibrate the toluene volume displacement before running the test.The (global) porosity was calculated from the measured bulk and particle densities as follows:

Properties of simulated porous bed and comparison with experimental data
The DEM model generated a porous bed consisting of 18,188 particles with sizes randomly distributed between 5.5 mm and 7.5 mm (Fig. 5).Based on the measured mass and volume of the porous bed, as well as the particle density and size, the total number of particles in the test bin was estimated 18,498.In other words, the total number of particles in the simulation was 1.7 % less than the measured value.It should be noted that estimation of the total number of particles in the bin was based on the assumption that all kernels were spherical and had an average diameter of 6.5 mm (the average of 5.5 and 7.5 mm used in simulations).However, the measured geometric mean diameter of particles was 6.6 mm (as discussed in section 2.1).Using the measured diameter of 6.6 mm, the total number of particles in the bin was estimated to be 17,670, which was 2.9 % less than the simulated value.
The simulated global (overall) porosity before vibration was 0.401, while the measured value was 0.398 ± 0.006 (± standard deviation, n = 6).The simulated porosity was within the 95 % CI (Confidence Interval) of the measured porosity (0.391 to 0.404).The simulated minimum porosity was 0.379, whereas the measured value was 0.366 ± 0.001 after vibration (at 1.0g intensity).The relative change (decrease) in porosity after vibration was 5.5 % in simulations, which was lower than that in measurements (8.0 %).Liu C. (2008) reported relative changes between 0.7 % to 6.1 % in porosity of soybean in a model bin after vibration.The simulated local porosity, calculated from the PFC 3D built-in function in a small region within the porous bed, varied with the location.Specifically, the local porosity before vibration was 0.402, 0.395, 0.395, 0.393, and 0.389 at the five selected locations (B1-B5), respectively.The average of the five locations was close to the measured value (0.395 vs. 0.398).Local porosity also varied with both vibration duration and intensity, which is discussed in the next sections.

Variation of porosity during vibration
Fig. 6 shows typical curves of simulated local porosity during vibration.The minimum porosity was reached at about 9 s for all locations.The local fluctuations in porosity were apparent, and they were attributed to the cyclic nature of vibration which compressed and relaxed the porous material during each vibration cycle (An X.Z. et al., 2005).The magnitudes of these local fluctuations could be as high as 0.018 (e.g., from 0.378 to 0.396 in porosity).The overall trend was that local porosity decreased as vibration progressed, but not monotonically (large global fluctuations were apparent).Taking location B2 as an example, the porosity increased quickly from the initial level of 0.395 to 0.430 in 0.04 s, and then dropped to 0.370 at about 0.3 s (Fig. 6).This observed pattern was similar  2005) observed that macro-pores (pores large enough to accommodate a particle) existed in loosely packed sphere assemblies.Therefore, high porosity could exist locally where a macro-pore exists.During vibration, particles re-arrange themselves and macro-pores would be destructed (An X.Z. et al., 2005).However, it was also possible that some macro-pores were created at an instant through particle re-arrangement during vibration, resulting in high local porosity.These instant macro-pores might not be stable and would be destructed quickly by vibration.During vibration, destruction and creation of macro-pores might occur repeatedly, resulting in fluctuations in local porosity; however, more macro-pores would be destructed than created, leading to a general trend of decrease in local porosity, as well as gradual decrease in the global porosity.A similar pattern of fluctuations in local porosity were observed by Liu C. (2008) in his experiments to investigate the vibration effect on porosity of soybean.He used a high speed camera to record movements of soybean kernels in a transparent bin made of Plexiglas.Based on the recorded video images, he tracked three particles that formed a triangle and calculated the area of this triangle at different times during vibration.He considered that changes in the triangle area reflected changes in local porosity (i.e., increasing triangle area reflected increasing in porosity, and vice versa).The relative change in local porosity was then estimated as S/S 0 , where S 0 was the initial area of triangle and S the triangle area at a given time.It was apparent that the relative change in local porosity fluctuated during vibration, with an overall  trend of decrease, similar to what was observed in his study (Fig. 7).

Effects of vibration intensity on porosity
Vibration intensity had noticeable effect on porosity during vibration (Fig. 8).The higher the intensity, the greater the fluctuations.Among the five intensity levels simulated (0.45g, 0.91g, 1.81g, 2.72g, and 3.62g), 1.81g resulted in the lowest porosity (Fig. 9).In other words, there was a "critical" vibration intensity (1.81g) that resulted in a maximum reduction in porosity, and lower or higher vibration intensities generated less dense porous structure.A similar observation was reported by An X. Z. et al. (2009).They measured changes in packing density of an assembly of spheres with diameter of 5.02 mm in a 229.70-mm diameter container when subjected to vibration between 0.25g and 5.0g, and observed that the highest packing density occurred in the mid-range vibration intensities between 2.0g to 2.5g.Theoretically speaking, if the acceleration of particles in a porous bed is greater than 1g in the vertical direction, the particles could become "afloat" once in every vibration cycle, which could result in a high probability of creating a large void into which another particle can slide, thus providing the best opportunity to form a dense structure.Therefore, lower vibration intensity (< 1g) would produce less reduction in porosity.However, it should be noted that a porous bed is a discrete system and not all particles vibrate at the same frequency and amplitude as the excitation applied on the containing structures (the bin bottom in this study).Energy transfer in the porous bed when subjected to vibration is extremely complicated, however it is a reasonable assumption that not all vibration energy from the vibration source applied at the bin bottom would be transferred to the particles in the porous bed because of energy dissipation by collision and friction among particles and between particle and the containing structures.Energy dissipation in granular materials during vibration has been shown to be significant in many studies in the area of particle impact damping (Friend R.D. and Kinra V.K., 2000).Energy dissipation offers a plausible conjecture for the reason why the maximum reduction in porosity occurred at 1.8g instead of 1g.This also means that the critical vibration intensity that produces the maximum reduction in porosity would vary with many factors that influence energy transfer in porous beds, including particle size and shape, particle density, friction between particles, friction between particles and containing structures, restitution coefficient, and even porosity itself, which dictates the distances (clearances) between particles and between particles and containing structure.
An X.Z.et al. ( 2009) explained that too high energy input would over-excite particles in porous beds, which would lead to less dense structures.While low amplitude vibration induces subtle changes in the microstructure of a granular material, resulting in higher bulk density, strong vibration leads to fluid-like behavior of the granular material, resulting in lower bulk density (Rosato A.D. et al., 2002). Hunt M.L. et al. (1994) conducted experiments to study the behavior of a granular material subjected to vertical vibration at acceleration levels from 1g to 5.5g.They observed that the volume of the granular material began to expand at acceleration of 2g.
8 Simulated average local porosity during vibration at intensities of 0.45g, 1.81g, and 3.62g.This meant that vibration energy at low intensity was not sufficient to cause significant particle re-arrangement, resulting in little change in porosity.Liu C. et al. (2010) used a high speed camera to study vibration of individual particles in a bin subjected to vertical vibration and observed that all particles in the bin vibrated as a continuum with the same frequency and amplitude as the excitation at low intensities (< 1g), whereas individual particles vibrated "irregularly" (not following the harmonic excitation) when the vibration intensities were high (> 1g).

Effect of vibration on airflow path
Vibration changed the shape and length of airflow path (Fig. 10).Qualitatively speaking, the deviation of flow path (the exit point) from the entry point was less before vibration than after vibration.This was probably because higher flow resistance in more dense porous beds after vibration made it more difficult for flow to go upwards directly.Nwaizu C. and Zhang Q. (2015) used imaging techniques (with smokes) to study airflow through grain bulks subjected to vibration and observed that reduction in porosity after vibration caused more lateral "branching" of flow when air entered the grain bulks vertically.In other words, more air was forced to flow in non-vertical directions after vibration.It was interesting to note that while there was less degree of deviation of flow path from the entry point before vibration, but more sharp local turns was observed.This became apparent when comparing the flow paths between 0g and 1.81g vibration intensities.Specifically, many small sharp local turns could be seen in the 0g path, whereas long sections of relatively smooth line existed in the 1.81g path (Fig. 10).
Both the deviation of flow path from the entry point and the local turns would result in longer flow paths (higher tortuosity).To quantify the differences in airflow path among different vibration intensities, the simulated tortuosity values were compared (Fig. 11).The tortuosity decreased after vibration generally, but no clear pattern could be identified.A maximum reduction of 9 % was observed for vibration intensity of 1.81g and a minimum 2 % for 2.72g.Many studies have shown that tortuousness of flow through porous media decreases with porosity (Pisani L., 2011).As discussed earlier, the lowest porosity occurred at 1.81g.Therefore, it was expected that flow path for 1.81g vibration would be less tortuous.This also meant that sharp local turns (dominant before vibration) had more effect on the tortuousness of flow path than did the deviation (dominant after vibration) from the entry point.Although the maximum reduction in tortuosity corresponded to the minimum porosity at 1.81g, the pattern of variation in tortuosity with vibration intensity did not quite follow that for porosity shown in Fig. 9.
The pore throat width (the equivalent diameter) of the widest local flow path was calculated for each tetrahedron unit associated with a global flow path (Fig. 12).It was observed that the throat width varied randomly between 1.3 and 4.6 mm, with one exception before vibration-a high value of 6.8 mm was observed.A pore throat of 6.8 mm was apparently greater than the average particle diameter.This meant that a macro-pore, defined as a pore large enough to accommodate a particle (An X.Z. et al., 2005), existed in the porous bed before vibration, but no macro-pores were observed after vibration at any of the  five vibration intensities studied.
To further examine the effect of vibration intensity on pore throat, the distributions of throat width were plotted (Fig. 13).The highest percentage of pore throats was in the 3-4 mm range before vibration, while the highest percentage after vibration was in the 2-3 mm range.The percentage of pore throats in the rages of 1-2 mm and 2-3 mm before vibration was much lower than that after vibration, whereas the percentage in 4-5 mm before vibration was much higher than that after vibration.These observations indicated that vibration reduced the pore throat width.On average, the pore throat width was 3.3, 2.6, 2.8, 2.7, 2.7, and 2.7 mm for 0g (before vibration), 0.45g, 0.91g, 1.81g, 2.72g, and 3.62g, respectively.While the average pore throat before vibration was greater (22 %) than those after vibration, the differences in pore throat width among the five different vibration intensities were negligible (< 0.2 mm).

Conclusions
(1) The discrete element model was capable of predicting some critical pore structure parameters for flow through porous beds, including porosity, tortuosity, pore throat width, and flow path geometry (length and shape).
(2) The influence of vibration on porosity was clearly demonstrated.A vibration intensity of 1.81g resulted in the lowest porosity (6 % lower than that before vibration).Lower vibration intensity did not have enough energy to cause much change in porosity, while higher intensity over-excited particles, producing less dense packing than did 1.81g.(3) Local porosity fluctuated markedly during vibration, probably due to repeated formation and destruction of instant (unstable) macro-pores (defined as pores large enough to accommodate a particle).(4) Tortuousness of flow path was noticeably affected by vibration.Tortuosity of the porous bed before vibra-tion was higher (2 %-9 %) than that after vibration.(5) Pore throats larger than the average particle diameter (due to macro-pores) existed in the porous bed before vibration, but were not observed after vibration.(6) Vibration reduced the pore throat width.The highest percentage of pore throats was in the range of 3-4 mm before vibration, and reduced to 2-3 mm after vibration.The reduction in average pore throat width was 22 % after vibration.

Fig. 1 A
Fig. 1 A Schematic of simulated bin showing the locations (B1-B5) where local porosity was calculated.

Fig. 2
Fig. 2 A tetrahedron as the basic building block in a porous bed.(a) four particles associated with a tetrahedron unit, (b) open space (a pore throat) in a tetrahedron unit.

Fig. 3
Fig. 3 Replacement of sharp turn angle between two local flow paths with arc.Fig. 4 Model bin test system (the bin was 0.28 m in height and filled to a depth of 0.25 m).
to the experimental observations by Nowak E.R. et al. (1998), but did not fully agree with the DEM simulation result reported by An X.Z.et al. (2005).They reported that packing fraction (opposite to porosity) increased monotonically with time of vibration (local fluctuations also existed).It should be noted that it was not clear if the porosity reported by An X.Z.et al. (2005) was the local or global porosity.Given that the pore structure in a porous bed varies from location to location, one would expect that local porosity varies with location.An X.Z.et al. (

Fig. 6
Fig. 6 Simulated local porosity during vibration at locations B1 and B2, and the average porosity of all five locations for vibration intensity of 1.8g.

Fig. 7
Fig. 7 Relative change in area of triangle formed by three particles during vibration at 25 Hz frequency and 0.45 mm amplitude (based on data reported by(Liu C., 2008).

Fig. 9
Fig. 9 Effect of vibration intensity on porosity.

Fig. 10
Fig. 10 Typical simulated airflow paths (as a series of associated tetrahedron units) for different vibration intensities.Fig. 11 Variation of simulated tortuosity with vibration intensities (each data point represents the average of three simulated paths and the error bar indicates the minimum and maximum values).

Fig. 12
Fig. 12 Pore throat width of local flow paths from the bin bottom (tetrahedron no. 1) to the top under different vibration conditions.