# Proper Orthogonal Decomposition Analysis and Dispersion Characteristics of Resonant Acoustic Flow

Authors: Xiaopeng Wang; Shifu Zhu; Song Chen; Ning Ma; Zhe Zhang
Journal: Shock and Vibration
 (2020)

---

## Abstract

The investigation on the flow field and mixing characteristics of resonant sound mixing is of great significance for the dispersion mixing of superfine materials. In order to simulate the flow field and dispersion characteristics of resonant acoustic mixing, a gas-liquid-solid three-phase flow model based on the coupled level-set and volume-of-fluid (CLSVOF) and discrete particle model (DPM) was established. The CLSVOF model solves the gas-liquid interface, and the DPM model tracks the particle position. Then, the particle image velocimetry (PIV) experiment was performed using a self-made resonance acoustic hybrid prototype under different oscillation accelerations, and the radial velocity distribution between the experiment and simulation was compared. Finally, the proper orthogonal decomposition (POD) is used to decompose the flow field under different oscillation accelerations and fill levels, and the energy distribution law and the energy structure of different scales are extracted. The results show that the energy of the instantaneous flow field of the resonant sound is mainly concentrated in the low-order mode, and a close relationship was revealed between the energy distribution law and dispersion behavior of particles. The larger the small-scale coherent structures distribute, the more energy it has and the more favorable it is for fast and uniform dispersion.


---

## Body

## 1. Introduction

As a new and efficient mixing technology, resonant acoustic mixing (RAM) has been applied in many fields such as nanomaterials, additive manufacturing, and chemical reactions. RAM generates violent fluid motion by seeking and operating at the “resonant condition” of the mechanical system at all times. Like all physical objects, a RAM system has a specific resonant frequency. This frequency around 60 Hz is constant when at rest but is ever-changing when in operation. RAM devices impart unique and differing physical effects on the mixing ingredients; research studies on the practical applications and the mixing mechanism have become highlights.American researchers conducted a RAM experiment of PBX explosives and propellants at the China Lake Military Base and found that the mixed materials of the process had higher uniformity than the mixing of conventional extruders [1]. Osorio et al. [2] observed changes in material properties of the final blend and examined variations in particle size, powder flow properties, and hydrophobicity of the blend. Their studies have shown that acceleration, blending time, and total energy input had a significant impact on powder blend and tablet properties. Vandenberg and Wille [3] optimized the mixing efficiency of RAM using ultrahigh performance concrete (UHPC) through acceleration curve profiles and workability spread flow tests and compared RAM with a top paddle mixer with regard to its effects on workability and flow, as well as compression and bending strength properties. Ende et al. [4] used RAM to mix the active pharmaceutical ingredient (API) compound and conformer at 60 Hz frequency to induce conversion to cocrystals without grinding media and produced CBZ : NCT cocrystals successfully. Park et al. [5] compared RAM technology with traditional ball-milling technology, and the results show that RAM technology can produce high-performance SOFC components efficiently and ultrafast. It can be found that the above studies are all based on experiments, and some research studies focus on numerical simulation of the flow field in vertically vibrated liquid column. Michalchuk et al. [6] studied the RAM-induced cocrystallisation process by using synchrotron X-ray powder diffraction for the first time and the phase diagram of the reaction between nicotinamide and carbamazepine in the presence of a small amount of water, and their studies provide a new perspective for understanding the role of various experimental parameters in traditional mechanical chemistry using liquid-assisted grinding techniques. Despite much research as above, there is still a lack of knowledge of the mixing mechanisms, the effects of operation parameters on dispersion process. With the development of computer technique, great progress has been made on studying fluid motion using the numerical simulation method. Buchanan et al. [7] developed a new procedure for bringing gas-liquid-solid three phases into contact to promote exchange phenomena and derived an equation which permits prediction of the minimum frequency for cyclic migration for the air-water system accurately. Friesen et al. [8] used level set and volume-of-fluid interface tracking algorithms to predict the experimentally determined bubble translation behavior accurately and confirmed the linear dependence of the bubble translation amplitude on the container amplitude. Xu et al. [9] investigated the effects of solid particles with very low holdup on the bubbling dynamics in a gas-liquid bubble column with single gas nozzle with coupled DPM and VOF. In order to clarify the influence of the combination of liquid meniscus and multiple droplets on the movement of liquid meniscus, the coupled level set and volume-of-fluid method (CLSVOF) and the continuum surface force (CSF) model were used by Wang et al. Then, they investigated the dynamic behavior of liquid meniscus movement, which is consistent with the experimental results [10].Numerical simulation method widely used in the flow field will generate a huge number of data, and the POD method has precisely unique advantages when dealing with massive data. The basic feature of the POD method is that a certain number of eigenmodes can be extracted from a series of instantaneous flow fields, and the flow field information is expressed by the eigenmode series expansion method. Hence, POD technique can be applied to understand the characteristics of complex flow and the relationship between flow field characteristics and mixing effects deeply. In the field of turbulence research, Lumley first introduced this method and used it to identify coherent structures in the flow field [11]. In the direct POD method, the spatial correlation matrix is difficult to solve this problem directly when the number of spatial points is large. Sirovich proposed the new mathematical processing method of the snapshot POD method, which is generally applicable to cases where the number of space points is much larger than the number of samples [12–15]. Owing to its reduced order and feature extraction characteristics, POD methods are increasingly being used to study complex flow phenomena. For example, Xiao et al. [16] obtained numerical simulations of molten iron flow fields in three different desulfurization processes and used POD to analyze flow field data. The analysis concluded that eccentric agitation is most conducive to desulfurization. Lehwald et al. [17] used particle image velocimetry (PIV) to measure the flow field of a static mixer and used POD to study the flow field structure and flow field visualization that are most closely related to the mixed field. Wang et al. [18] used POD to study the effect of fluid viscoelasticity on the turbulence characteristics of a double-oscillation grid. Wen et al. [19] used TR-PIV to measure the flow field of the submerged jet and free-surface interaction. Then, they used POD to decompose the measured flow field and explained the interaction between the vortex and the free surface.Oscillation acceleration and fill level have great influence on the dispersion efficiency in the process of dispersion, and there is a close relation between the flow field in the mixing vessel and them; hence, it makes sense to study the flow field in the mixing vessel. In this article, the information of the flow field with different conditions was extracted by POD technique from a series of transient velocity data that were obtained by numerical simulation. Then, the characteristics of POD modes and energy spectrum were compared and analyzed. The relationship between the energy distribution law and dispersion behavior of particles is finally elucidated.


## 2. Mathematical Model and Simulation Methods

The free surface and the motion of the gas-liquid phase both have strong influence on the entire flow field in the dispersion process; therefore, mathematical models included in this CFD simulation are the CLSVOF method accounting for the gas-liquid phase interface and the DPM model dealing with the dispersion process of micron particles [20, 21]. The coupling between the CLSVOF model and the DPM model was realized by the momentum exchange between the continuous phase and the particles.


### 2.1. CLSVOF Model

The gas-liquid phase flow in the mixing container contains physical phenomena such as shearing, collision, and crushing, and the gas-liquid phase interface is complicated. The VOF [22] approach is an interface-tracking method of multiple incompatible phases based on a fixed Euler grid, and it uses phase volume fractions to characterize the sharing of different fluids in each computational grid unit. Different fluid components share a set of momentum equations throughout the computational domain.Based on the VOF model, the CLSVOF [23] model used in this paper calculates the volume fraction through the level-set function. The CLSVOF model makes up the inconsistency of the VOF model reconstruction process and the inaccurate orientation of the surface. At the same time, it corrects the gas-liquid phase interface and the mixing flow field. The level-set function is defined as a signed distance to the interface. Accordingly, the interface is the zero level-set, σx,t, and can be expressed as Γ=xσx,t=0 in a two-phase flow system: (1)σx,t=+dI,ifx∈the first phase,0,ifx∈Γ,−dI,ifx∈the secondary phase,(2)∂φ∂t+∇⋅uφ=0,where dI is the distance from the interface.The evolution of the level-set function Ф can be given in a similar form as to the VOF model:(3)∂ϕ∂t+∇⋅uϕ=0.Also, the momentum equation can be written as(4)∂ρu∂t+∇⋅ρuu=−∇p+∇⋅μ∇u+∇uT−Fsf+ρg,where u is the fluid phase velocity and Fsf is the force arising from surface tension effects.In view of the lack of variable conservation and the equality of sharp interface in the calculation process of the level-set method, the CLSVOF model corrects level-set function by volume fraction.


### 2.2. Discrete Particle Model

The position of a discrete phase particle was calculated by integrating the force balance on the particle without mutual interaction, which was enforceable in a Lagrange coordinate system [24]. Particle force differential equation represents a balance between the particle inertia and external forces acting on the particle and can be written (for the x direction in Cartesian coordinates) as(5)dupdt=FDu−up+gxρp−ρρp+Fx,where up is the particle velocity, gx is the component of the gravitational acceleration in the x direction, ρ is the fluid density, ρp is the particle density, Fx is an additional acceleration term in the particle force balance that can be important under special circumstances, mainly including pressure gradient force Fpr=upρ/ρp∂u/∂x, virtual mass force Fvm=1/2ρ/ρpd/dtu−up, and Saffman’s lift force (Fsl=2Kv1/2ρdij/ρpdpdlkdklu−up). The virtual mass and pressure gradient forces are important unless the density of the fluid is much lower than the density of the particles as is the case for liquid/solid particles in gaseous flows ρ/ρp≪1, [25] and Saffman’s lift force is intended for small particle Reynolds numbers and is recommended only for submicron particles. Hence, these forces were taken into consideration in this paper. FDu−up is the drag force per unit particle mass, and(6)FD=18μρpdp2⋅CDRe24,where μ is the molecular fill level of the fluid and dp is the particle diameter. Re is the Reynolds number, which is defined as(7)Re=ρdpup−uμ.For spherical particles, the drag coefficientCD is(8)CD=a1+a2Re+a3Re2,where a1, a2, and a3 are constants over several ranges of Re given by Morsi and Alexander [26].


### 2.3. Boundary Conditions and Numerical Simulation Strategies

The geometry of gas-liquid-solid flow simulation system is a cylinder. The fluid domain of the half section of the container axis is taken for simulation owing to the axisymmetric container. The no-slip boundary condition (BC) was imposed on the container walls, and the container axis was considered as axisymmetric BC. The particles define the details in a file and are added to the liquid surface. More detail for these simulations can be found in Table1.Table 1
The parameters of the numerical model.
Scheme
Parameter
Value
Structure
Container diameter (mm)
84
Container height (mm)
94
Material
Liquid density (kg/m3)
998.2
Liquid dynamic viscosity (Pa·s)
0.001
Gas density (kg/m3)
1.225
Gas dynamic viscosity (Pa·s)
1.79 × 10−5
Particle density (kg/m3)
1500
Particle diameter (μm)
15
Gas-liquid tension coefficient (N/m)
0.073
Operating
Vibration acceleration (g)
20–40
Fill level (%)
70–90
Vibration frequency (Hz)
50
Particle number
1000
g (m/s [2])
9.8Considering turbulent eddies in the dispersion process and near-wall treatment, the renormalization group (RNG)k-ε turbulent model was selected. Enhanced wall treatment was adopted for near-wall treatment, and Y+ was close to 1 [27]. The geometry was meshed with GAMBIT 2.4.6. All zones were meshed with structured grids with the maximum grid size of 1 mm, and the area near wall was refined with 16 elements considering the use of wall functions. This generated grids of 5,800 cells for the whole computation domain. The material object and grids are shown in Figure 1. Such a meshing method undergoes the check of grid independence.Figure 1
Real products, diagram of the vessel, and meshes of computational domain, respectively.The numerical simulation was carried out using the ANSYS Fluent software. The sliding grid boundary condition was used to simulate the vertical vibration of the container, and the Fluent UDF macro DEFINE_TRANSIENT_PROFILE (V, time) was used to define the periodic motion. The pressure-based Navier–Stokes solution algorithm was used to separate and solve model equations, pressure-velocity coupling was the pressure-implicit with the splitting of operator (PISO) method, the pressure discretization method was pressure staggering option (PRESTO!), and the momentum discretization method was second-order upwind scheme. The Geo-Reconstruct discretization scheme was used to solve volume fraction equations for the VOF explicit scheme. In order to prevent the solution failure caused by the large global courant number, the time step size was 0.0001 s. Computer’s configuration is as follows: CPU: Xeon E5-2699 32 cores, memory: 128G, and network environment: Windows 7 64 bit. Eight cores were used for parallel computation. Numerical simulation result is shown in Figure 2. Due to the combination of surface tension and wall adhesion, the liquid phase shows a tendency of “sticking to the wall.” The liquid splits into droplets of different sizes. Bubbles are drawn into the liquid and then collide, coalesce, and break. The consistent distribution of particles and bubbles confirms the important role of gas-liquid distribution in dispersion.Figure 2
Numerical simulation result att = 0.1 s, 0.3 s, and 0.5 s, respectively. (a) The flow pattern. (b) The position of the particles.
(a)
(b)


### 2.4. Numerical Simulation Verification

The experiments were carried out based on the PIV system as shown in Figure3. The experimental platform was mainly composed of the RAM prototype, a double-pulsed laser, high-resolution CCD camera, laser sheet optics, timing circuit, power supply, image-sampling computer, and image processing software Dynamic Studio (version 3.14). The camera was operated at 2048 × 2048 pixels with a framing rate of 15 Hz. Each interrogation area included 32 × 32 pixels (with 50% overlap in each direction). Postprocessing work was performed with Tecplot 360. The velocity vector obtained by the PIV experiment is shown in Figure 4, which included the radial velocity u and the axial velocity v. Since the axial velocity was affected by the vibration velocity of the container itself, the radial velocity was selected to compare the experiment with simulation. In order to avoid the influence of turbulence randomness on the result, the radial velocity was averaged after taking the absolute value:(9)u¯i=1n∑j=1nuij,where i, j is the axial and radial sequence number of the PIV measuring point, n is the number of measuring points of the same axial height, and uij is the radial velocity at the measuring point.Figure 3
Schematic representation of the resonant acoustic mixer (left) and the PIV system (right).Figure 4
Velocity vector on the axial section of 20 g, 30 g, and 40 g, respectively.Figure5 shows the comparison of radial velocity between the experiment and simulation at different accelerations. The radial velocity increases first then decreases with the increase of axial height. Finally, it increases and reaches its maximum at a height of about 60 mm. The container is divided into three regions a, b, and c according to the trend of the radial velocity. In region a, turbulence near the liquid surface decreases in intensity as it propagates downward, so the radial velocity increases with the increase of axial height. In region b, the influence of the fluctuation of the liquid surface is continuously reduced as the height increases, so the radial velocity decreases with the increase of axial height. Region c is closer to the wall of the vessel, and the radial velocity increases with the increase of axial height owing to the stronger impact between the liquid phase and the wall. In addition, as the acceleration increases, the radial velocity at the same axial height generally shows an increasing trend. This is mainly due to the liquid level fluctuation, the liquid phase and the wall impact and the bubble generation are more intense with the increase of the acceleration, which increases the turbulence intensity. Owing to the limitations of processing and installation accuracy, the movement of the container is not strictly axial, which leads to error in the experimental and simulation results. In general, the simulation results have the same trend with the experimental results, and the numerical simulation is credible.Figure 5
Simulations and experiments of radial velocity in gas-liquid-solid flows.


## 2.1. CLSVOF Model

The gas-liquid phase flow in the mixing container contains physical phenomena such as shearing, collision, and crushing, and the gas-liquid phase interface is complicated. The VOF [22] approach is an interface-tracking method of multiple incompatible phases based on a fixed Euler grid, and it uses phase volume fractions to characterize the sharing of different fluids in each computational grid unit. Different fluid components share a set of momentum equations throughout the computational domain.Based on the VOF model, the CLSVOF [23] model used in this paper calculates the volume fraction through the level-set function. The CLSVOF model makes up the inconsistency of the VOF model reconstruction process and the inaccurate orientation of the surface. At the same time, it corrects the gas-liquid phase interface and the mixing flow field. The level-set function is defined as a signed distance to the interface. Accordingly, the interface is the zero level-set, σx,t, and can be expressed as Γ=xσx,t=0 in a two-phase flow system: (1)σx,t=+dI,ifx∈the first phase,0,ifx∈Γ,−dI,ifx∈the secondary phase,(2)∂φ∂t+∇⋅uφ=0,where dI is the distance from the interface.The evolution of the level-set function Ф can be given in a similar form as to the VOF model:(3)∂ϕ∂t+∇⋅uϕ=0.Also, the momentum equation can be written as(4)∂ρu∂t+∇⋅ρuu=−∇p+∇⋅μ∇u+∇uT−Fsf+ρg,where u is the fluid phase velocity and Fsf is the force arising from surface tension effects.In view of the lack of variable conservation and the equality of sharp interface in the calculation process of the level-set method, the CLSVOF model corrects level-set function by volume fraction.


## 2.2. Discrete Particle Model

The position of a discrete phase particle was calculated by integrating the force balance on the particle without mutual interaction, which was enforceable in a Lagrange coordinate system [24]. Particle force differential equation represents a balance between the particle inertia and external forces acting on the particle and can be written (for the x direction in Cartesian coordinates) as(5)dupdt=FDu−up+gxρp−ρρp+Fx,where up is the particle velocity, gx is the component of the gravitational acceleration in the x direction, ρ is the fluid density, ρp is the particle density, Fx is an additional acceleration term in the particle force balance that can be important under special circumstances, mainly including pressure gradient force Fpr=upρ/ρp∂u/∂x, virtual mass force Fvm=1/2ρ/ρpd/dtu−up, and Saffman’s lift force (Fsl=2Kv1/2ρdij/ρpdpdlkdklu−up). The virtual mass and pressure gradient forces are important unless the density of the fluid is much lower than the density of the particles as is the case for liquid/solid particles in gaseous flows ρ/ρp≪1, [25] and Saffman’s lift force is intended for small particle Reynolds numbers and is recommended only for submicron particles. Hence, these forces were taken into consideration in this paper. FDu−up is the drag force per unit particle mass, and(6)FD=18μρpdp2⋅CDRe24,where μ is the molecular fill level of the fluid and dp is the particle diameter. Re is the Reynolds number, which is defined as(7)Re=ρdpup−uμ.For spherical particles, the drag coefficientCD is(8)CD=a1+a2Re+a3Re2,where a1, a2, and a3 are constants over several ranges of Re given by Morsi and Alexander [26].


## 2.3. Boundary Conditions and Numerical Simulation Strategies

The geometry of gas-liquid-solid flow simulation system is a cylinder. The fluid domain of the half section of the container axis is taken for simulation owing to the axisymmetric container. The no-slip boundary condition (BC) was imposed on the container walls, and the container axis was considered as axisymmetric BC. The particles define the details in a file and are added to the liquid surface. More detail for these simulations can be found in Table1.Table 1
The parameters of the numerical model.
Scheme
Parameter
Value
Structure
Container diameter (mm)
84
Container height (mm)
94
Material
Liquid density (kg/m3)
998.2
Liquid dynamic viscosity (Pa·s)
0.001
Gas density (kg/m3)
1.225
Gas dynamic viscosity (Pa·s)
1.79 × 10−5
Particle density (kg/m3)
1500
Particle diameter (μm)
15
Gas-liquid tension coefficient (N/m)
0.073
Operating
Vibration acceleration (g)
20–40
Fill level (%)
70–90
Vibration frequency (Hz)
50
Particle number
1000
g (m/s [2])
9.8Considering turbulent eddies in the dispersion process and near-wall treatment, the renormalization group (RNG)k-ε turbulent model was selected. Enhanced wall treatment was adopted for near-wall treatment, and Y+ was close to 1 [27]. The geometry was meshed with GAMBIT 2.4.6. All zones were meshed with structured grids with the maximum grid size of 1 mm, and the area near wall was refined with 16 elements considering the use of wall functions. This generated grids of 5,800 cells for the whole computation domain. The material object and grids are shown in Figure 1. Such a meshing method undergoes the check of grid independence.Figure 1
Real products, diagram of the vessel, and meshes of computational domain, respectively.The numerical simulation was carried out using the ANSYS Fluent software. The sliding grid boundary condition was used to simulate the vertical vibration of the container, and the Fluent UDF macro DEFINE_TRANSIENT_PROFILE (V, time) was used to define the periodic motion. The pressure-based Navier–Stokes solution algorithm was used to separate and solve model equations, pressure-velocity coupling was the pressure-implicit with the splitting of operator (PISO) method, the pressure discretization method was pressure staggering option (PRESTO!), and the momentum discretization method was second-order upwind scheme. The Geo-Reconstruct discretization scheme was used to solve volume fraction equations for the VOF explicit scheme. In order to prevent the solution failure caused by the large global courant number, the time step size was 0.0001 s. Computer’s configuration is as follows: CPU: Xeon E5-2699 32 cores, memory: 128G, and network environment: Windows 7 64 bit. Eight cores were used for parallel computation. Numerical simulation result is shown in Figure 2. Due to the combination of surface tension and wall adhesion, the liquid phase shows a tendency of “sticking to the wall.” The liquid splits into droplets of different sizes. Bubbles are drawn into the liquid and then collide, coalesce, and break. The consistent distribution of particles and bubbles confirms the important role of gas-liquid distribution in dispersion.Figure 2
Numerical simulation result att = 0.1 s, 0.3 s, and 0.5 s, respectively. (a) The flow pattern. (b) The position of the particles.
(a)
(b)


## 2.4. Numerical Simulation Verification

The experiments were carried out based on the PIV system as shown in Figure3. The experimental platform was mainly composed of the RAM prototype, a double-pulsed laser, high-resolution CCD camera, laser sheet optics, timing circuit, power supply, image-sampling computer, and image processing software Dynamic Studio (version 3.14). The camera was operated at 2048 × 2048 pixels with a framing rate of 15 Hz. Each interrogation area included 32 × 32 pixels (with 50% overlap in each direction). Postprocessing work was performed with Tecplot 360. The velocity vector obtained by the PIV experiment is shown in Figure 4, which included the radial velocity u and the axial velocity v. Since the axial velocity was affected by the vibration velocity of the container itself, the radial velocity was selected to compare the experiment with simulation. In order to avoid the influence of turbulence randomness on the result, the radial velocity was averaged after taking the absolute value:(9)u¯i=1n∑j=1nuij,where i, j is the axial and radial sequence number of the PIV measuring point, n is the number of measuring points of the same axial height, and uij is the radial velocity at the measuring point.Figure 3
Schematic representation of the resonant acoustic mixer (left) and the PIV system (right).Figure 4
Velocity vector on the axial section of 20 g, 30 g, and 40 g, respectively.Figure5 shows the comparison of radial velocity between the experiment and simulation at different accelerations. The radial velocity increases first then decreases with the increase of axial height. Finally, it increases and reaches its maximum at a height of about 60 mm. The container is divided into three regions a, b, and c according to the trend of the radial velocity. In region a, turbulence near the liquid surface decreases in intensity as it propagates downward, so the radial velocity increases with the increase of axial height. In region b, the influence of the fluctuation of the liquid surface is continuously reduced as the height increases, so the radial velocity decreases with the increase of axial height. Region c is closer to the wall of the vessel, and the radial velocity increases with the increase of axial height owing to the stronger impact between the liquid phase and the wall. In addition, as the acceleration increases, the radial velocity at the same axial height generally shows an increasing trend. This is mainly due to the liquid level fluctuation, the liquid phase and the wall impact and the bubble generation are more intense with the increase of the acceleration, which increases the turbulence intensity. Owing to the limitations of processing and installation accuracy, the movement of the container is not strictly axial, which leads to error in the experimental and simulation results. In general, the simulation results have the same trend with the experimental results, and the numerical simulation is credible.Figure 5
Simulations and experiments of radial velocity in gas-liquid-solid flows.


## 3. Proper Orthogonal Decomposition

Following the snapshot POD ideas of Sirovich [12], the steps for calculating the POD modes and the corresponding coefficients can be derived as follows:(10)∫ΩRx,x′φx′dx′=λφx,where Ω represents the spatial domain x=x,y∈Ω considered and λ and φ are the coefficients and eigenmodes, respectively. The kernel R(x, x’) is the two-point correlation function defined as(11)Rx,x′=1N∑n=1Nux,tnuTx′,tn,where u(x, tn) is the two-dimensional data of velocity fields and N is the number of snapshots.Assuming that the eigenmodes can be written in terms of combination of the original dataset as(12)φx=∑n=1NAtnux,tn,and inserting (12) and (13) into (11), after some simplifications, it yields(13)CA=λA.MatrixC and eigenvector A in equation (14) are determined as follows:(14)C=Ci,j=1NuTx,tiux,tii,j=1,…,N,(15)A=Atn.By eigen decomposition of matrixC in equation (14), a set of eigenvalues and eigenvectors is solved and then used to compute the POD modes according to equation (13). The POD modes are obtained by standard normalization of the eigenvectors as(16)φix=∑n=1NAitnux,tn∑n=1NAitnux,tni=1,…,N.In this particle, the method was applied on 2D experimental datasets, which correspond to N planes consisting ofk × m = M data points. All velocity vectors can be arranged in matrix U as follows:(17)U=u1u2…uN=u11u12…u1N⋮⋮⋮⋮uM1uM2…uMNv11v12…v1N⋮⋮⋮⋮vM1vM2…vMN.Meanwhile, matrixC is formed as(18)C=UTU.Arrange the eigenvalues obtained by equation (14) in the descending order, and we obtain(19)λ1≥λ2≥⋯≥λN≥0.Thenth eigenvalue λn represents the average turbulent kinetic energy of the nth eigenmode, and the energy ratio En of the nth eigenmode is as follows:(20)En=λn∑i=1Nλi.In addition, the cumulative percentage to the total turbulent kinetic energy TKEn of the first n eigenmodes is as follows:(21)TKEn=∑i=1nEi.


## 4. Results and Discussion

As it is mentioned in relevant literature, the sample space required for POD analysis has an impact on the analysis structure. Based on the snapshot number independence verification, the number of snapshots of all the later analysis conditions in this paper isN = 100.


### 4.1. POD Energy Spectrum

The POD energy spectrum is plotted in Figure6, which shows the energy distribution En of the nth mode and the cumulative contribution TKEn of the first n modes for different fill levels at oscillation acceleration = 40 g. It can be clearly seen that, for the 90% fill level, energy is less distributed in high-order modes, but for the 70% or 80% fill level, energy is distributed widely over a large number of POD modes. It is indicated that fill level has a significant impact on flow field characteristics. Different order POD modes mean different scale coherent structures; therefore, there is a corresponding relationship between the POD energy spectrum and the flow complexity. The broad POD energy spectrum represents that the flow field needs to be described by more coherent structures. The POD energy spectrum of the 70% fill level is the broadest; the 80% fill level flow takes the second place, and the POD energy spectrum of the 90% fill level is the narrowest. Meanwhile, we should note that the difference between the 80% fill level flow and 70% fill level flow is far less than the difference between the 80% fill level flow and 90% fill level flow. It means that the complexity of flow decreases with increasing fill level and changes slowly when the fill level is large enough. Similarly, Figure 7 shows the energy distribution En of the nth mode and the cumulative contribution TKEn of the first n modes at oscillation acceleration = 20 g, 30 g, and 40 g under 70% fill level. It is obvious that the energy gradually migrates to higher-order POD modes as the oscillation acceleration increases, and the complexity of flow has a positive correlation with the oscillation acceleration.Figure 6
Energy distribution of POD modes under different fill levels at 40 g oscillation acceleration.Figure 7
Energy distribution of POD modes at different oscillation accelerations under 70% fill level.According to the POD principle, the eigenvalue is given the meaning of the turbulent kinetic energy, and thenth eigenvalue represents the energy of the nth POD mode. So, the energy distribution En of the nth mode and the cumulative contribution TKEn of the first n modes for the 70%, 80%, and 90% fill level at oscillation acceleration = 20 g, 30 g, and 40 g should be further analyzed, listed in Table 2. From Table 2, it can be seen that the first POD mode captures 19.25%, 15.00%, 8.84%, 9.13%, and 9.93% of turbulent kinetic energy in the five cases respectively, and for the 63th mode, which are 0.23%, 0.30%, 0.42%, 0.39%, and 0.06%, respectively. It suggests that overlarge fill level has a significant inhibitory effect on high-order POD modes, and oscillation acceleration is more closely related to low-order POD modes. From Table 2, it is found that 40, 54, 63, 61, and 46 POD modes are required in the five cases, respectively. The number of POD modes that capture more than 90% of the energy can reflect the concentration of turbulent kinetic energy in the flow field. This is mainly owing to low oscillation acceleration of the mixing machine cannot provide enough energy to dissipate turbulence to a small scale; hence, the energy is concentrated in the low-order modes, and the excessive fill level reduces the volume of the gas in the container and the gas holdup of the liquid. At the same time, the space for the free movement of the liquid is also reduced, and the impact between the liquid and wall is weakened, which will hinder the generation of small-scale vortex structures represented by high-order modes. Based on the above comparison, the flow with different oscillation accelerations or fill levels exhibits a huge difference.Table 2
The energy distributionEn of the nth mode and the cumulative contribution TKEn of the first n modes.
n
70% fill levela = 20 g
70% fill levela = 30 g
70% fill levela = 40 g
80% fill levela = 40 g
90% fill levela = 40 g
E
n
TKEn
E
n
TKEn
E
n
TKEn
E
n
TKEn
E
n
TKEn
1
19.25
19.25
15.00
15.00
8.84
8.84
9.13
9.13
9.93
9.93
40
0.46
90.56
0.63
82.83
0.75
77.34
0.69
79.21
0.48
82.83
46
0.38
91.67
0.52
86.20
0.64
81.42
0.58
82.93
0.31
90.03
54
0.29
92.89
0.40
90.05
0.54
86.14
0.50
87.24
0.19
92.01
61
0.24
94.74
0.33
92.27
0.45
89.55
0.41
90.35
0.08
92.63
63
0.23
95.21
0.30
92.90
0.42
90.41
0.39
91.15
0.06
93.28In POD theory, the original flow field can be reconstructed by a series of POD modes, and these POD modes should capture more than 90% of the energy [28, 29]. It indicates that the flow under 70% and 80% fill level at oscillation acceleration = 40 g has more freedom degrees than that of the other three cases.


### 4.2. POD Modes

In POD, a series of flow field data will be decomposed onto a number of time-invariant basis functions (POD modes), and different POD modes represent spatial structures of different scales. The high-order POD modes characterize the vortex structure with large scale and strong turbulent kinetic energy; on the contrary, the low-order POD modes characterize the vortex structure with small scale. We can easily analyze the characteristics of vortexes of various scales in the flow field through POD modes instead of direct analysis of the original flow field. POD modes are still composed of the velocity vector essentially, and we express the POD modes as the vorticity field by using Tecplot 360 so as to show the vortex structure more clearly. Figures8–12 show vorticity fields depicted by the 1st, 4th, 16th, and 64th POD mode for five different cases, and the red dotted line indicates the position of the liquid level.Figure 8
The single POD mode at 20 g oscillation acceleration under 70% fill level.Figure 9
The single POD mode at 30 g oscillation acceleration under 70% fill level.Figure 10
The single POD mode at 40 g oscillation acceleration under 70% fill level.Figure 11
The single POD mode at 40 g oscillation acceleration under 80% fill level.Figure 12
The single POD mode at 40 g oscillation acceleration under 90% fill level.Combined with the flow pattern, it can be clearly seen that the large-scale vortexes gradually dissipate into small-scale vortexes with the increase of the POD mode order. Liquid is not strong enough to split under the viscous and gravitational forces and collide directly with the wall of the container due to the low acceleration as shown in Figure8. Since the vortexes of the first POD mode appear at the liquid level, with the increase of the mode order, large-scale vortex structures gradually dissipate into small-scale vortex structures, which reflect the pulsation characteristics of the higher velocity field. The vortexes of the 64th POD mode mainly distribute above the liquid level, which will result in a slow transfer rate between the liquid phase regions. As the oscillating acceleration increases from 20 g to 30 g, liquid can break through the action of viscous force and gravity and directly collide with the wall surface of the container. With the impact of liquid, bubbles of various scales inside liquid go through a process of generation—development—change (deformation, polymerization, and fragmentation)—collapse (as shown in Figure 2), and turbulence pulsation intensifies; hence, the small-scale vortexes of the 64th POD mode are more widespread as shown in Figure 9, and vorticity intensity increases from 1.2 to 4.2. Increased turbulence pulsation will promote dispersion. When the oscillating acceleration reaches 40 g, the small-scale vortexes are evenly distributed throughout the container, and it reaches the optimal state (whole-field mixing state) of dispersion at this time.The actual production process usually uses high fill level to improve production efficiency. As the fill level increases from 70% to 80%, the free space in the container is compressed, and the process of bubbles from generation to collapse is restrained. There is no significant change in the distribution range of small-scale vortexes and vorticity intensity. It means that the increase of the fill level has little effect on the distribution characteristics of vortexes. While as the fill level comes to 90%, the inhibitory effect becomes apparent. The distribution range and intensity of the small-scale vortexes decrease sharply owing to the weakening of the gas-liquid impact strength. It suggests that there is an optimal fill level about 80% for dispersion.


### 4.3. Particle Disperse Process Analysis

In order to evaluate the uniformity of particle distribution obtained by the DPM model, a feature point distribution uniformity evaluation method in the field of image processing was employed. The calculation domain was divided into 3 × 3 subareas, the number of particles in each area was one sample data, and the number of particles in all areas of the container became a sample set. The standard deviation of the sample sets represented difference in the particle numbers of different areas. The smaller the standard deviation, the more uniform the particle distribution within the container. The calculation process is as follows:(22)σ=1n∑i=1nNi−N¯2,where n is the number of subareas, Ni is the number of particles in the ith region, and N− is the average number of particles in n regions.Figure13 shows the effects of oscillating acceleration and fill level on particle dispersion efficiency. It can be seen that the case under 70% fill level at 40 g has the highest dispersion efficiency; the case under 80% fill level at 40 g takes the second place, whereas the gap between the 70% and 80% fill level is very small; then the case under 70% fill level at 30 g and the case under 90% fill level at 40 g; the case under 70% fill level at 20 g disperses the slowest. Considering that the fill level is as large as possible to improve the production efficiency, the most favorable condition for particle dispersion is 80% fill level and 40 g. Combined with the above POD mode analysis, it is not difficult to see that the more energy high-order POD modes contain, the faster particles disperse. That is to say, the intensity and distribution range of the small-scale vortexes are closely related to the dispersion process.Figure 13
The standard deviation of particle distribution under five different cases.


## 4.1. POD Energy Spectrum

The POD energy spectrum is plotted in Figure6, which shows the energy distribution En of the nth mode and the cumulative contribution TKEn of the first n modes for different fill levels at oscillation acceleration = 40 g. It can be clearly seen that, for the 90% fill level, energy is less distributed in high-order modes, but for the 70% or 80% fill level, energy is distributed widely over a large number of POD modes. It is indicated that fill level has a significant impact on flow field characteristics. Different order POD modes mean different scale coherent structures; therefore, there is a corresponding relationship between the POD energy spectrum and the flow complexity. The broad POD energy spectrum represents that the flow field needs to be described by more coherent structures. The POD energy spectrum of the 70% fill level is the broadest; the 80% fill level flow takes the second place, and the POD energy spectrum of the 90% fill level is the narrowest. Meanwhile, we should note that the difference between the 80% fill level flow and 70% fill level flow is far less than the difference between the 80% fill level flow and 90% fill level flow. It means that the complexity of flow decreases with increasing fill level and changes slowly when the fill level is large enough. Similarly, Figure 7 shows the energy distribution En of the nth mode and the cumulative contribution TKEn of the first n modes at oscillation acceleration = 20 g, 30 g, and 40 g under 70% fill level. It is obvious that the energy gradually migrates to higher-order POD modes as the oscillation acceleration increases, and the complexity of flow has a positive correlation with the oscillation acceleration.Figure 6
Energy distribution of POD modes under different fill levels at 40 g oscillation acceleration.Figure 7
Energy distribution of POD modes at different oscillation accelerations under 70% fill level.According to the POD principle, the eigenvalue is given the meaning of the turbulent kinetic energy, and thenth eigenvalue represents the energy of the nth POD mode. So, the energy distribution En of the nth mode and the cumulative contribution TKEn of the first n modes for the 70%, 80%, and 90% fill level at oscillation acceleration = 20 g, 30 g, and 40 g should be further analyzed, listed in Table 2. From Table 2, it can be seen that the first POD mode captures 19.25%, 15.00%, 8.84%, 9.13%, and 9.93% of turbulent kinetic energy in the five cases respectively, and for the 63th mode, which are 0.23%, 0.30%, 0.42%, 0.39%, and 0.06%, respectively. It suggests that overlarge fill level has a significant inhibitory effect on high-order POD modes, and oscillation acceleration is more closely related to low-order POD modes. From Table 2, it is found that 40, 54, 63, 61, and 46 POD modes are required in the five cases, respectively. The number of POD modes that capture more than 90% of the energy can reflect the concentration of turbulent kinetic energy in the flow field. This is mainly owing to low oscillation acceleration of the mixing machine cannot provide enough energy to dissipate turbulence to a small scale; hence, the energy is concentrated in the low-order modes, and the excessive fill level reduces the volume of the gas in the container and the gas holdup of the liquid. At the same time, the space for the free movement of the liquid is also reduced, and the impact between the liquid and wall is weakened, which will hinder the generation of small-scale vortex structures represented by high-order modes. Based on the above comparison, the flow with different oscillation accelerations or fill levels exhibits a huge difference.Table 2
The energy distributionEn of the nth mode and the cumulative contribution TKEn of the first n modes.
n
70% fill levela = 20 g
70% fill levela = 30 g
70% fill levela = 40 g
80% fill levela = 40 g
90% fill levela = 40 g
E
n
TKEn
E
n
TKEn
E
n
TKEn
E
n
TKEn
E
n
TKEn
1
19.25
19.25
15.00
15.00
8.84
8.84
9.13
9.13
9.93
9.93
40
0.46
90.56
0.63
82.83
0.75
77.34
0.69
79.21
0.48
82.83
46
0.38
91.67
0.52
86.20
0.64
81.42
0.58
82.93
0.31
90.03
54
0.29
92.89
0.40
90.05
0.54
86.14
0.50
87.24
0.19
92.01
61
0.24
94.74
0.33
92.27
0.45
89.55
0.41
90.35
0.08
92.63
63
0.23
95.21
0.30
92.90
0.42
90.41
0.39
91.15
0.06
93.28In POD theory, the original flow field can be reconstructed by a series of POD modes, and these POD modes should capture more than 90% of the energy [28, 29]. It indicates that the flow under 70% and 80% fill level at oscillation acceleration = 40 g has more freedom degrees than that of the other three cases.


## 4.2. POD Modes

In POD, a series of flow field data will be decomposed onto a number of time-invariant basis functions (POD modes), and different POD modes represent spatial structures of different scales. The high-order POD modes characterize the vortex structure with large scale and strong turbulent kinetic energy; on the contrary, the low-order POD modes characterize the vortex structure with small scale. We can easily analyze the characteristics of vortexes of various scales in the flow field through POD modes instead of direct analysis of the original flow field. POD modes are still composed of the velocity vector essentially, and we express the POD modes as the vorticity field by using Tecplot 360 so as to show the vortex structure more clearly. Figures8–12 show vorticity fields depicted by the 1st, 4th, 16th, and 64th POD mode for five different cases, and the red dotted line indicates the position of the liquid level.Figure 8
The single POD mode at 20 g oscillation acceleration under 70% fill level.Figure 9
The single POD mode at 30 g oscillation acceleration under 70% fill level.Figure 10
The single POD mode at 40 g oscillation acceleration under 70% fill level.Figure 11
The single POD mode at 40 g oscillation acceleration under 80% fill level.Figure 12
The single POD mode at 40 g oscillation acceleration under 90% fill level.Combined with the flow pattern, it can be clearly seen that the large-scale vortexes gradually dissipate into small-scale vortexes with the increase of the POD mode order. Liquid is not strong enough to split under the viscous and gravitational forces and collide directly with the wall of the container due to the low acceleration as shown in Figure8. Since the vortexes of the first POD mode appear at the liquid level, with the increase of the mode order, large-scale vortex structures gradually dissipate into small-scale vortex structures, which reflect the pulsation characteristics of the higher velocity field. The vortexes of the 64th POD mode mainly distribute above the liquid level, which will result in a slow transfer rate between the liquid phase regions. As the oscillating acceleration increases from 20 g to 30 g, liquid can break through the action of viscous force and gravity and directly collide with the wall surface of the container. With the impact of liquid, bubbles of various scales inside liquid go through a process of generation—development—change (deformation, polymerization, and fragmentation)—collapse (as shown in Figure 2), and turbulence pulsation intensifies; hence, the small-scale vortexes of the 64th POD mode are more widespread as shown in Figure 9, and vorticity intensity increases from 1.2 to 4.2. Increased turbulence pulsation will promote dispersion. When the oscillating acceleration reaches 40 g, the small-scale vortexes are evenly distributed throughout the container, and it reaches the optimal state (whole-field mixing state) of dispersion at this time.The actual production process usually uses high fill level to improve production efficiency. As the fill level increases from 70% to 80%, the free space in the container is compressed, and the process of bubbles from generation to collapse is restrained. There is no significant change in the distribution range of small-scale vortexes and vorticity intensity. It means that the increase of the fill level has little effect on the distribution characteristics of vortexes. While as the fill level comes to 90%, the inhibitory effect becomes apparent. The distribution range and intensity of the small-scale vortexes decrease sharply owing to the weakening of the gas-liquid impact strength. It suggests that there is an optimal fill level about 80% for dispersion.


## 4.3. Particle Disperse Process Analysis

In order to evaluate the uniformity of particle distribution obtained by the DPM model, a feature point distribution uniformity evaluation method in the field of image processing was employed. The calculation domain was divided into 3 × 3 subareas, the number of particles in each area was one sample data, and the number of particles in all areas of the container became a sample set. The standard deviation of the sample sets represented difference in the particle numbers of different areas. The smaller the standard deviation, the more uniform the particle distribution within the container. The calculation process is as follows:(22)σ=1n∑i=1nNi−N¯2,where n is the number of subareas, Ni is the number of particles in the ith region, and N− is the average number of particles in n regions.Figure13 shows the effects of oscillating acceleration and fill level on particle dispersion efficiency. It can be seen that the case under 70% fill level at 40 g has the highest dispersion efficiency; the case under 80% fill level at 40 g takes the second place, whereas the gap between the 70% and 80% fill level is very small; then the case under 70% fill level at 30 g and the case under 90% fill level at 40 g; the case under 70% fill level at 20 g disperses the slowest. Considering that the fill level is as large as possible to improve the production efficiency, the most favorable condition for particle dispersion is 80% fill level and 40 g. Combined with the above POD mode analysis, it is not difficult to see that the more energy high-order POD modes contain, the faster particles disperse. That is to say, the intensity and distribution range of the small-scale vortexes are closely related to the dispersion process.Figure 13
The standard deviation of particle distribution under five different cases.


## 5. Conclusion

In this article, the numerical simulation was carried out with different oscillation accelerations and fill levels, and the simulation data were verified by PIV experiment. Then, the simulation data were processed by POD technique, and we compared POD characteristics and particle dispersion efficiency under different cases. The main conclusions are drawn as follows.Through the analysis of the POD energy spectrum, the POD energy spectrum of the 70% fill level at 40 g is the broadest. It indicates that the complexity of flow decreases with increasing fill level and decreasing oscillating acceleration and changes slowly when the fill level is large enough.Compared POD modes with particle dispersion efficiency in five different cases, the small-scale vortexes represented by high-order POD modes have close relationship with particle dispersion efficiency, that is, the more energy high-order POD modes contain, the faster particles disperse.Based on the above analysis of the flow field, POD, and particle dispersion process, we reveal the effect of the small-scale vortexes on the dispersion process and obtain the optimal parameters for particle dispersion.



---

*Source: 5068042-2020-03-04.xml*