Non-equilibrium effect modeling and particle velocity estimation for particulate flow in porous media

The complex dynamics of fluid and particles flowing through pore space demands some relaxation time for particles to catch up with fluid velocity, which manifest themselves as non-equilibrium (NE) effects. Previous studies have shown that NE effects in particle transport can have significant consequences when relaxation time is comparable to the characteristic time associated with the fluid flow field. However, the existing models are lacking to account for this complicated relation between particles and fluid. In this paper, we adapt the general form of the harmonic oscillation equation to describe NE effects in the particulate flow system. The NE effect is evaluated by solving coupled mass balance equations with computational fluid dynamic (CFD) techniques using COMSOL Multiphysics®. A simplified straight-tube model, periodic converging–diverging tube model, and SEM image of a real pore network are applied in NE analyses. The results indicate that the time variation of the NE effect complies with the theory of stability. Two key parameters of the oscillator equation are amplitude (A) and damping ratio (ζ), where the former represents the magnitude of NE and the latter is an indication of flow path geometry. The velocity equations for particle transport in different flow path geometries are derived from the proposed NE equation, offering a quick estimation of particle velocity in the particulate flow. By conducting simulation on the SEM image of a real pore structure, the equivalent radii of the pores where particles move through were obtained. The outcome of this work can shed light on explaining the complex NE effects in porous media. The generalized equation to model NE can help temporarily decouple particle transport equations from fluid equations, facilitating much advanced particulate flow modeling problems in the large-scale problems.


Introduction
Particle transport behavior is the issue that petroleum engineers attempted to solve for many years because particulate flow systems exist in a wide range of petroleum-related processes, for example, the injection of seawater during water flooding, invasion of filtrate during drilling, and micro-enhanced oil recovery (Yuan et al., 2016). The formation damage caused by fine-particle retention and detachment highly affects the productivity and injectivity of wells (Bedrikovetsky et al., 2011). A large amount of fine production may also result in equipment erosion, flowline plugging, and other potential hazards (Marquez et al., 2014). In the asphaltene deposition problem, asphaltene particles can change to accumulate or adsorb onto the pore surfaces, which leads to pore plugging in the reservoir and finally affects the flow rate within the wellbore (Davudov et al., 2018). During the proppant placement process, the performance of proppant packs in hydraulic fracture plays a significant role in fracture conductivity and production behavior (Fan et al., 2018). It is believed that the uniform distribution of proppants provides the highest fracture conductivity. However, it is very challenging to sustain uniform distribution because of proppant settlement and embedment (Yu and Sepehrnoori, 2013). Therefore, it is extremely important for petroleum engineers to understand particle transport behavior in order to apply it in the aforementioned petroleum processes.
Engineers have derived many mathematical models and conducted several experiments to investigate the particle transport mechanisms in particulate systems. In the previous research studies of particle transport in particulate flow through porous media, it has been discovered that the particle velocity does not necessarily have the same velocity as its carrier fluid when moving through pore networks. It had been observed that the velocity of particles near the pore surface is less than the velocity of their carrier fluid (Yuan and Shapiro, 2010). In a vertical domain, particles tend to accelerate in a short period before reaching the terminal velocity (Awad et al., 2021). Moreover, the particle settling velocity is found to be dropped in an oscillatory flow (Amaratunga et al., 2021). Researchers conducted a nanoparticle flooding experiment and observed that nanoparticles may temporarily stick in the core (Hu et al., 2016). The reason is that when nanofluids flow from pores to throats, the flowing channel area becomes narrow, causing an increase in the nanofluid velocity. The fluid velocity will become faster than nanoparticles. It causes nanoparticles to accumulate at the entrance of the pore throat (Sun et al., 2017). When addressing the issue of the large-scale proppant transportation in fractures, petroleum engineers rely on two methods: (1) the velocity model saying that viscosity produces the influence of the proppant load on the fluid and (2) the mixed phase model in which the proppant and fluid are considered as different phases with different velocities (Roy et al., 2015). The first model is not sophisticated because it is only applicable for high viscosity fluids, and it cannot describe the correct particle behavior in low viscosity fluids. The second model (mixed phase model) tracks particles and fluids separately; therefore, it cannot capture the complex behavior between them.
In this study, the complex dynamic between the particle and fluid in porous media and how it influences particle transport behavior will be investigated. The particle velocity equation will be derived from the proposed NE equation, offering a new approach to decouple the particle equation from the fluid equation without solving the complicated multiphase model. The processes are organized as follows: Section 2 introduces the method to solve coupled mass balance equations with CFD techniques. Moreover, the NE equation is proposed. Sections 3 and 4 present the modeling of NE in different flow path geometries. The particle velocity equations are decoupled from fluid equations by incorporating the proposed NE equation. In Section 5, the application of the NE equation in actual porous media is addressed. Finally, the conclusion is drawn.

Solving particulate flow using the CFD technique
There are two computational fluid dynamic (CFD) approaches that treat particles and fluids separately when dealing with the transport problem in hydraulic fracture or porous media. The first one is called the Lagrangian approach, also referred to as the particle tracking method or the discrete phase model (DPM) (Kong et al., 2016;Song and Park, 2019;Wang et al., 2019). This method treats the fluid as a continuum phase by solving the Navier-Stokes (N-S) equation and then tracks each particle as a discrete phase by coupling them with the flow field. The limitation of this method is that the volume of the disperse phase cannot be employed because the DPM assumes particle loading is low compared with the whole domain, and the particle-particle interaction is not being considered. However, if the DPM is coupled with discrete element method (DEM), the particle-particle and particle-wall interactions can be captured. DEM is a method of applying Newton's law to particles to obtain their motions. The second method is the Euler-Euler approach, which treats each phase as an interpenetrating continuous phase. It solves a set of momentum and mass balance equations for both solid and liquid phases. The solid phase is considered as the continuum fluid phase of the particle. The Euler-Euler model is Basic schematic of straight tube generated using COMSOL.

Frontiers in Earth Science
frontiersin.org the most complex and computationally intensive among the multiphase models. In this study, the particulate flow was solved using the particle tracking method using COMSOL Multiphysics ® . Particles or droplets were considered as rigid particles. The interaction between the particle and fluid was obtained by coupling Newton's law with the flow field. The governing equation for particle transport in the fluid flow is where ρ is the fluid density, u is the fluid velocity, μ is fluid viscosity and F is the additional force acting on the fluid. In this model, it is assumed that all particles are naturally buoyant; therefore, the gravity and buoyancy forces can be ignored, and the drag force (F D ) is the dominant factor in determining the particle trajectories. The governing equation for particle motion is where v is the particle velocity and m P is the mass of the particle. The general form of drag force is where τ P is the relaxation time. Based on the Stokes drag law, the particle relaxation time is defined as In this study, the NE effect is defined as a function of particle velocity (v p ) and local fluid velocity (v f ), which can be written as

Adapting the harmonic oscillator equation to describe NE
In the mechanical vibration system, the harmonic oscillator equation with damping is driven by solving Newton's second law equation, which is where A is called the oscillation amplitude and ζ is called the damping ratio. The system behavior depends on the value of the damping ratio (ζ).
• Overdamped (ζ >1): The system returns to steady state without oscillating • Critically damped (ζ =1): The system returns to the steady state as quickly as possible without oscillating • Underdamped (ζ <1): The system oscillates with the amplitude gradually decreasing to zero The underdamped condition is what we used for the NE effect modeling. The general solution for ζ < 1 can be written as This study is driven by the given hypothesis: the linear theory of stability can explain non-equilibrium evolution in particulate systems through the general form of the harmonic oscillation equation. The NE parameter is proposed to be Combining Eq. 5 with Eq. 8, the adapted form of the harmonic oscillator equation yields to be

FIGURE 2 (A)
Velocity profiles in the cross section and (B) particle trajectories in the straight-tube model. Color scale represents the magnitude of fluid and particle velocities.

Model setup
The basic geometry schematic of the straight-tube model is shown in Figure 1. A straight tube with a uniform radius of 100 μm and a total length of 250 μm was modeled. In order to reduce the computational cost, the geometry was cut to a quarter of the original tube along the symmetry line.
The fluid was given a velocity of 0.004 m/s at the inlet boundary. The outlet boundary was set at p=0, referring to the inlet boundary. Fifty particles were injected from different positions at the inlet boundary to the outlet boundary. The details of the simulation model are listed in Table 1.

Fluid velocity profile and particle trajectory
The velocity profile and particle trajectory are illustrated in Figure 2. The flow velocity distribution across the tube follows the character of Poiseuille flow. With the only effect of the drag force, particles move along the flow streamline.
Using the Hagen-Poiseuille Law, the fluid velocity at any given radius (r) inside the tube can be expressed as where v max is the maximum fluid velocity inside the tube, and R is the radius of the tube. In this case, v max is around 0.006 m/s, and R is equal to 100 μm (0.0001m). Plugging them into Eq. 10, the velocity distribution inside the tube is v f r ( ) 0.006 1 − r 1 × 10 −4 2 . (11)

Non-equilibrium parameter determination
The NE parameter (1 − vp vf ) was evaluated using COMSOL. Simulation results were matched with the harmonic oscillator equation using MATLAB ® to obtain the magnitude (A) and damping ratio (ζ) values. The P50 value on cumulative distribution functions (CDF) of A and ζ were recorded separately to represent the most likely value.
The NE parameter as a function of time for one of 50 particles is plotted in Figure 3A. At an early stage (when particle is just injected into the tube), there is velocity difference between the particle and fluid. In that time period, particle velocity is always less than fluid velocity. After a certain time period, the particle and fluid reach the same velocity, which can be described as an equilibrium state. The adapted harmonic oscillation equation (Eq. 9) was implemented to match simulation results using MATLAB. The curve-fitting result for a single particle is shown in Figure 3B.
The same curve-fitting process was performed for all 50 particles. The histograms and cumulative distribution functions for two key factors of NE, A and ζ, are presented in Figure 4. The P50 values of A and ζ at their cumulative distribution functions are 0.025 and 0.977, respectively.

Decouple particle equation from fluid equation
The P50 values of ω 0 and φ on their cumulative distribution functions were determined to be 305 and 100, respectively. Plugging the P50 values of all coefficients into Eq. 9, the following equation is obtained: Combining Eq. 11 with Eq. 12, the particle velocity for any given radius (r) and time (t) for the case of μ 1 cp is determined to be FIGURE 3 (A) NE parameter as a function of time solved using COMSOL, and (B) curve fitting for the NE effect of a single particle using MATLAB.

Model setup and simulation results
In a real situation, the complex roughness geometries of fractures or porous media are often characterized by triangular, rectangular, trapezoidal, bell, sinusoidal, and random shapes (Dejam et al., 2018). In this work, for simplicity, a periodical converging-diverging-shaped tube is used to characterize porous media. The basic schematic of the converging-diverging-shaped tube model is demonstrated in Figure 5A. Two diverging-shaped tubes and one convergingshaped tube are connected with each other. Each tube has a maximum radius of 100 μm and a minimum radius of 50 μm, as well as a tube length of 250 μm. Fifty particles were injected at the same velocity with the initial fluid velocity from the inlet boundary to the outlet boundary. The details of the simulation are same as those of the straight-tube model except that the

Fluid velocity profile
The velocity profile v(r, z) in a tube whose cross-section area is gradually changing with the distance along the flow direction remains parabolic (Bahrami et al., 2008). The schematic velocity profile is shown in Figure 6. Therefore, the axial velocity needs to be described as where u m (z) is the average velocity at the axial location z and a(z) is the tube radius at the location z. By applying the conservation of mass, we obtain where a 0 is the tube radius at the initial point. Combining Eq. 14 and 15, the axial velocity equation is For a convergent-shaped tube, where R min , R max , and L are the minimum tube radius, maximum tube radius, and tube length, respectively. Similarly, for a divergentshaped tube, By combining Eq. 16 with Eq. 17, we get the generalized velocity equation in a converging tube, which is By combining Eq. 16 with Eq. 18, we get the generalized velocity equation in a diverging tube, which is Plugging R min 5 × 10 −5 m, R max 1 × 10 −4 m, L 2.5 × 10 −4 m, and 2u m,0 =0.0045 m/s into Eq. 19, the fluid equation for any given radius (r) and axial length (z) inside the convergent is v f r, z ( ) 4.5 × 10 −3 1 × 10 −4 1 × 10 −4 − z (22)

FIGURE 6
Schematic of varying the cross-section tube.
Frontiers in Earth Science frontiersin.org

Non-equilibrium parameter determination
Following the same steps as the straight-tube model, the NE parameters were evaluated using COMSOL. Simulation results were matched with the harmonic oscillator equation to obtain the magnitude (A) and damping ratio (ζ) values. The P50 values on cumulative distribution functions (CDF) of A and ζ were recorded separately.
The NE parameter as a function of time is shown in Figure 7A According to the figure, the particle experiences three main stages depending on the flow path geometry. The flow geometry for the first stage is a divergent-shaped flow path. The NE parameter is less than zero, which means that particle velocity exceeds fluid velocity when the tube radius is getting bigger. At an early time, there is a big difference between the particle velocity and fluid velocity. At late time, particle velocity gets close to fluid velocity, indicating they are close to the equilibrium state. After the first stage, the flow path geometry switches to a convergent shape. The NE parameter is larger than zero. It implies that particle velocity is always less than fluid velocity. In contrary to system behavior in the divergent flow path, the particle and fluid velocity will never getting close to each other, which indicates that the non-equilibrium will always happen. The reason for the NE parameter drops at late time is because the flow path geometry transits from converting to a diverging pattern.
The harmonic oscillation equation was implemented to match simulation results for different stages separately. The initial time for all three stages was normalized to zero for curve fitting. The curve fitting result for a single particle in stage 2 is shown in Figure 7B. The A and ζ values were determined to be 2.4 × 10 −5 and -0.978. In frictional vibration, negative damping causes system instability (Chen, 2014). In this situation, it indicates that the NE parameter increases as a function of time. Since particle flowed in the converging tube, it can be concluded that the negative ζ value implies particle to move in convergent flow path geometry. Stage 1 and stage 3 had same particle-fluid behavior, and the fitting curve is shown in Figure 7C. The A and ζ values were determined to be 0.026 and 0.94. In this case, ζ equals to a positive number. It means that the NE parameter decreases over time. Since particle flowed in the diverging flow path, it can be concluded that the positive ζ value is an indication of particle moving in divergent flow pattern.
After matching all 50 particles' data, the histograms and cumulative distribution function for A and ζ in divergent and convergent flow patterns were obtained separately. The P50 values of A and ζ for the convergent flow pattern are determined to be 1.9 × 10 −5 and -0.894, respectively. As for divergent flow path geometry, P50 values for A and ζ are 0.015 and 0.926, respectively. The result complies with the single particle case where ζ is negative when the flow pattern FIGURE 7 (A) NE parameters as a function of time in the converging-diverging tube solved using COMSOL, and the curve fitting of NE for single particle injection for (B) stage 2, and (C) stages 1 and 3 using MATLAB.

Frontiers in Earth Science
frontiersin.org is convergent, whereas ζ is positive when the flow pattern is divergent shape. Moreover, the result indicates that the A value for the convergent flow pattern is smaller than that of the divergent flow pattern.

Decouple particle equation from the fluid equation
The particle transport equation is obtained by combining the fluid equation with the NE equation. All coefficients in the harmonic oscillation equation that are used to describe NE were obtained from P50 values on their cumulative distribution functions. The P50 values of A, ζ, ω 0 , and φ for convergent flow path geometry are 1.9×10 -5 , -0.894, 421, and 100, respectively.
Plugging P50 values of those coefficients into Eq. 9, the NE equation for the divergent flow pattern is obtained as By combining Eq. 22 with Eq. 25, the particle velocity equation for any given radius (r), axial location (z), and time (t) in the convergent tube is determined to be v p r, z, t ( ) 0.014 5 × 10 −5 5 × 10 −5 + z 5 2 1 − r 5 × 10 −5 + z In order to validate the accuracy of the particle velocity equation that decoupled from the fluid equation, the predicted particle velocities with different initial locations for convergent and divergent flow path geometries from the analytical particle velocity equation are compared with the simulation results solved using COMSOL Multiphysics ® . The results are shown in Figure 8. It is shown in the figure that there are decent matches between the analytical result and simulation result, which proves the accuracy of the model. 5 Non-equilibrium effects in the actual pore network

Model setup and simulation results
The model geometry that we used is a SEM image of a real pore network from the COMSOL application library. The binary image and imported geometry are shown in Figure 9A. The size of the image is 640 μm by 320 μm. The black region represents the pore space. The white region represents the rock grain. The single particle was injected from the inlet (right side) to outlet (left side) boundary in the same velocity as the fluid velocity. The particle and fluid properties are set to be the same as those of the straight-tube model. The fluid velocity distribution and single particle trajectory inside the pore structure are demonstrated in Figures 9B, C.

Pore size distribution from SEM images
The pore size distribution is determined using the MATLAB code designed by Rabbani et al. (2014). It is based on the watershed segmentation algorithm, which can detect and Assuming that each segment is a circle, the pore radius can be determined. The cut structure and pore size distribution are shown in Figure 10.

Non-equilibrium parameter determination
The NE effect for the single particle was evaluated using COMSOL. The NE as a function of time for a single particle injection is plotted in Figure 11A. Due to the complexity of the pore structure, the trends of NE variation over time are not identical. As discussed earlier, NE larger than zero indicates that particles move through the converging flow pattern, whereas NE less than zero implies that particles move in the diverging flow pattern. In addition, different-shaped pores are supposed to have different oscillation behaviors, which manifest as distinct curve shapes in the figure. Therefore, based on the flow pattern criteria and the shape of the NE curve, the NE curve can be divided into several stages. For example, from t=0 to t=0.01s, the NE is above zero, which is an indication of converging flow path geometry. Meanwhile, there are two distinct curve shapes, which means that the particle moves in two different pores.
After detecting the flow path and the curve shape, 20 stages were identified. Again, the adapted harmonic oscillation equation was  Frontiers in Earth Science frontiersin.org implemented to match the data for each stage. The initial time for each stage was normalized to zero for curve fitting. The oscillator amplitudes and damping ratios were obtained and are demonstrated in Table 2. The results present that in the converging flow pattern, ζ values are less than zero, which is consistent with the result obtained in the previous model. It means that the particle velocity is always less than the fluid velocity, and they can never achieve equilibrium with each other. In the diverging flow pattern, the result shows that ζ values are greater than zero, which is also consistent with the previous observation. As discussed earlier, particle velocity is greater than fluid velocity in the diverging flow path. Using all coefficients (A, ζ, ω 0 , and φ) obtained from the curve fitting, the matched NE curve for all stages can be obtained, and they are shown in Figure 11B. The fitted curves have good matches with the original NE curve, which are evaluated in simulation.

Equivalent tube radius determination for the pore structure
In order to have a better characterization of particle transport behavior in porous media, the motions of multiple particles in different size of pores need to be investigated. Simulations were performed for different numbers of particle injections. The simulation details are kept same as those of the straight-tube model.
The simulation was performed on 5-20 particle injections (Jin, 2018). There are some pore spaces that particles can never reach with only the effect of drag force. Moreover, it was observed that some particles share same paths inside the pore structure. The more particles are injected, the more repeated flow paths are observed. Therefore, some tube radii will be calculated more than once finally yielding an over-estimation of equivalent pore number. Therefore, in this case, five particles' trajectories have covered almost all possible flow paths that particles can reach. It provides a better characterization of the pore network, and the corresponding tube radius distribution is considered the actual pore spaces that particles move through in this pore network.
In order to determine the radius distribution of pores that particles move through, the MATLAB code was designed to identify different stages in which particles move through different sizes of pores, and then the proposed NE equation was applied to match the results for all stages to obtain the A values.
The rationale for defining an equivalent radius is that since each stage has its corresponding A value in multiple particles' injection  simulation, it is hypothesized that there is an equivalent tube radius that can give the same A value for each stage. In order to determine the equivalent pore radii for all stages, another converging-diverging model simulation, which has been discussed in section 4, was performed. According to the pore size distribution of the pore structure, the minimum and maximum pore radii are 4.81 μm and 43 μm, respectively. The converging-diverging models were set up with a fixed minimum tube radius of 5 μm, and the variable maximum tube radius ranged from 10 μm to 45 μm. The length of each tube was set at 50 μm. The proposed NE equation was applied again to match the simulation results for the convergent and divergent flow pattern separately to obtain A values in different flow patterns. The relationships between the maximum radius and A value for convergent and divergent flow patterns are shown in Figure 12.
The results indicate that there are approximately linear relationships between the tube radius and A value for different flow patterns.
Using the trend line equations, the maximum radius for each A value can be determined. The average value of maximum and minimum tube radii was approximated to be the equivalent tube radius of the corresponding stage. For the convergent flow pattern, the equivalent tube radius (r e ) is estimated to be r e R max + R min 2 79375.5A − 16.
For the divergent flow pattern, the equivalent tube radius (r e ) is estimated to be r e R max + R min 2 3116.7A − 2.3.
The estimated equivalent pore radius distribution for paths where five particles move through is shown in Figure 13.

Conclusion
In this study, the NE effect in the particulate flow system was investigated by developing particle transport models in the straight tube, periodic converging-diverging-shaped tube, and actual pore structure. This study was hypothesizing whether the linear theory of stability can explain NE evolution in particulate systems through the general form of the harmonic oscillation equation. Particle velocity equations for convergent and divergent flow patterns are decoupled

FIGURE 12
Corresponding A values for different maximum radii in the convergent flow pattern and divergent flow pattern.

FIGURE 13
Histograms for the estimated tube radius for the five-particle injection scheme.

Frontiers in Earth Science
frontiersin.org from fluid equations by incorporating the proposed NE equation. In this work, only particle-fluid interaction is taken into consideration; the drag force is the dominant factor, and the carrier fluid is a singlephase fluid. The main conclusions are drawn as follows: 1) The time variation of the NE effect complies with the theory of stability. The harmonic oscillation equation can be used to characterize the NE effect in the particulate flow system. 2) Two key parameters of the oscillator equation are amplitude (A) and damping ratio (ζ). The former represents the magnitude of NE and the latter is an indication of flow path geometry. 3) In divergent flow path geometry, the ζ value is between 0 and 1.
The NE effect decreases as a function of time. The ζ value is between 0 and -1 in convergent flow path geometry where the NE effect increases as a function of time implying that the particle velocity always remains less than the fluid velocity; hence, the system will never achieve an equilibrium state. 4) The derived particle equation offers a quick estimation of particle velocity for a given time and the location of porous media. 5) The flow simulations of the SEM image present consistent results with diverging and converging flow results. By conducting different numbers of particle injection schemes, it is found that five-particle injection offers a better estimation of actual pores that particles move through. The corresponding equivalent tube radii of complex pore geometries are obtained.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.