Kinetic modeling of fractal aggregate mobility

Multiphase flows composed of the gas carrier and particulate dispersed phases are ubiquitous in nature as well as in industrial applications. For example, these systems are found in cyclone separators used in separating granular particulates based on particle size differences,¹ in sandblasting used for surface smoothening,² and combustion reactors in power generation.³ In nature, these systems are everywhere from aerosol pollutants in the air to dust storms in deserts.⁴ Similar multiphase systems are also encountered in the specific operations of lunar landings⁵ where the lunar dust or regoliths are entrained in the plume generated by the lander^6,7 as well as in high-speed aircraft moving through dust storms or dusty air.⁸ The transport or mobility parameters used in modeling the momentum transfer between the phases are important in predicting the evolution of the multiphase system as a whole. These parameters describing the translational and rotational dynamics of the particulate phase are in turn functions of particulate geometry as well as gas phase parameters. Since the momentum transfer between the phases is a function of multiple variables, such as particle geometry, orientation of the flow vector and non-linear function of gas parameters at higher Reynolds numbers, theoretical or analytic determination of such coupling parameters becomes difficult.

The sphere is the most commonly studied particulate geometry. Extensive studies on this geometry at various flow conditions at different levels of particle Reynolds numbers and gas rarefaction (Knudsen number) are available by several authors.^9–17 The spherically symmetric geometry simplifies the mobility by neglecting the net momentum transfer in the flow lateral direction (i.e., lateral drag) and the induced moment (i.e., pitching torque). Another extensively considered particle geometry is the ellipsoidal or spheroidal geometry. Although the volume of work for this geometry is not as extensive as that of the spherical geometry, researchers, such as Taylor,¹⁸ Jeffery,¹⁹ Ouchene et al.,²⁰ Zastawny et al.,²¹ Livi et al.,²² and Chinnappan et al.,^23,24 have studied the mobility of ellipsoidal particles in continuum to free-molecular gas flow regimes. These regular geometries are asymmetric about certain orientations with respect to the incoming flow velocity vector, and hence, particulates experience lateral drag and pitching torque along these orientations. However, a large fraction of particulates in nature and industries have highly irregular shapes rather than being fully spherical or ellipsoidal.²⁵ These non-spherical, irregular, angular particulates are found in powders of crystalline materials, such as glass, sand, and coal powder.²⁶ A specific class of irregular particulates may be modeled by the aggregation of primary particles or monomers. These constituent monomers usually have regular geometric construct, such as spherical, ellipsoidal, or spheroidal. Since these aggregates have fractal-like geometry and are highly porous making, the analysis of their mobility parameters becomes further difficult.²⁷ 

Several researchers have studied the mobility of fractal aggregates (FA) using numerical, analytic, and experimental techniques. The mobility of these grains was primarily studied in the Stokes regime (⁠$Re≪1$⁠), where the viscous effects of the fluid or gas dominate over the inertial effects. Fillipov et al.²⁸ extended a multipole expansion technique developed in light scattering problems to determine the forces and moments generated on arbitrary aggregates of spherical monomers. Researchers, such as Bossis et al.²⁹ and Gastaldi et al.,³⁰ used Stokesian dynamics to determine these parameters numerically for FAs. Theoretical work by Saarloos³¹ used the Kirkwood–Risemann approximation to determine the hydrodynamic radius, an important transport parameter of FAs. The above-mentioned studies were performed for the low-speed gas continuum flow regime, i.e., in the Stokes regime where the flow governing equation is essentially linear in gas viscosity and speed, while the works of Chan et al.,³² Nakamura et al.,³³ and Mackowski et al.³⁴ studied the drag induced on aggregates of spherical primary particles in the slow-moving, free-molecular gas regime using the direct simulation Monte Carlo (DSMC) method for the gas and observed the increasing importance of pressure drag on net mobility with increasing gas rarefaction. More recently, Zhang et al.³⁵ used the same approach to determine the scalar friction factor for a broader range of Knudsen number or gas rarefactions covering near-continuum to free-molecular flows.

The effect of gas compressibility and rarefaction on the fractal aggregate mobility is an important aspect in areas, such as lunar landing where the soil grains with the structure of fractal aggregates are entrained in the lunar lander plume of a rarefied high-velocity gas. Conventional and contemporary numerical methods used to study the mobility parameters of such complex structures are continuum-based Navier–Stokes (NS) solvers that require mathematical modeling of the gas–surface interaction (slip and temperature jump) to emulate rarefaction at such micrometer length scales. DSMC has the advantage of being one of the most efficient high-fidelity physics-based methods, which resolve the slip conditions implicitly.

In this article, we present the general mobility characteristics of small, rigid, inert model fractal aggregates in high-speed gas flows (Mach number < 6) at particle Knudsen numbers corresponding to slip and transitional regimes. The effect of the orientation of the FA with the flow vector as well as the effect of gas rarefaction on particle mobility is studied in detail by employing our in-house high-fidelity gas kinetic DSMC code CHAOS. Due to the computational efficiency of the CHAOS code, the particulate mobility study will be useful in understanding and modeling highly compressible multiphase flows, such as dust-shock interactions through highly rarefied flows, such as planetary landing where the rocket plume interacts with the dusty planet surface. The remainder of the paper is organized as follows: a brief discussion on the DSMC method used to model the gas-particulate system and the associated simulation parameters are given Sec. II. First, the method is employed to study the gas-induced mobility of a single spherical monomer in Sec. III. The geometric description of the model FA used in this study and the mathematical parameters describing the same will be crucial in analyzing our simulation results provided in Sec. IV. Results from the detailed analyses of the mobility of complicated fractal structures are presented in Sec. V. Section VI compares the mobility parameters of a FA with an ellipsoidal approximation under the same conditions to understand the sensitivity of the inner porous construct of the aggregate on its mobility. The important research findings from these studies are presented in Sec. VIII.

The gas-particulate multiphase system is modeled with the gas-kinetic DSMC³⁶ method using the in-house, GPU-based code CHAOS.³⁷ The gas phase is modeled using coarse-grained DSMC gas-computational particles representing a large number of actual gas molecules. The evolution of these DSMC gas particles through streaming and collisions approximates the solution of the Boltzmann equation because all numerical parameters are specified. The streaming or motion of the DSMC particles is modeled as a deterministic process while the binary, inter-molecular collisions are modeled as a stochastic process. The decoupled nature of these processes facilitates flexibility in using higher simulation time step magnitudes for the method.³⁶ A detailed description of the DSMC algorithm and the implementation of the same in our in-house code CHAOS is given in Jambunathan et al..³⁷ The momentum exchange between the gas molecules and the particulate surface through molecular collisions generates fluid-induced forces and moments on the particulate. In this work, the gas–surface interactions are modeled on a molecular level through collisions of DSMC gas particles with the triangular surface elements of the particulate through the immersed boundary method (IBM).^37,38 This physics-based approach emulates the surface boundary conditions at a fundamental level compared to continuum-based methods³⁹ that require mathematical models for velocity-slip boundary condition. Note that the precise gas–surface interaction modeling ensures that we will predict the mobility parameters of particulate mobility as accurately as possible.

The interaction of gas molecules or DSMC simulated gas particles with the triangular particulate surface elements coupled to the gas phase through the immersed boundary method is modeled as diffuse reflections to emulate a particulate rough surface as shown in Fig. 1. In this approach, prospective surface collisional particles of velocity $v→i$ are detected using an efficient ray-tracing algorithm built into the CHAOS code. The particle is then reflected from the surface in a random vector direction with a velocity $v→r$ sampled from a half-Maxwellian velocity distribution corresponding to the surface temperature of the particulate panel.

The net momentum transfer $Δp→pan$ on a particulate surface panel due to $ncoll$ DSMC particle collisions per simulation time step is given as follows:

where m_g and FNUM represent the molecular mass of the gas and the number of real gas molecules approximating a single DSMC simulated particle, respectively.

The net steady-state force experienced by the aggregate due to the gas–surface momentum transfer, $F→net$⁠, is obtained by integrating the momentum transferred $Δp→pan$ along the surface as follows:

where N_pan is the total number of triangular panels representing the particulate surface. The net force in this equation is resolved along the three orthogonal axes (x, y, z), but later in the study, we will consolidate the three components into drag and lift forces along the flow direction and perpendicular to the flow direction, respectively.

Similarly, the net moment or the pitching torque, $M→net$⁠, generated on the FA is given as follows:

where $r→p,CM$ is the direction vector between the center of mass of the particulate and the centroid of the surface panel as indicated in Fig. 1. Note that in Eqs. (2) and (3), $Δt$ represents the DSMC time step and sampling is performed for N_sam time steps after the system reaches the steady-state to reduce the statistical noise associated with the DSMC method.

The particulates considered in this study are exposed to gas flow (in the positive Z direction) with different free-stream Mach numbers in the range (0.1 and 6.5) for two specific Knudsen numbers (Kn) 0.035 and 0.35 corresponding to near-continuum (slip) and transitional regimes, respectively. Note that the Knudsen number is the ratio of the gas mean-free path λ calculated using the hard-sphere model and the diameter of the spherical monomer $dp=2 μ$m.³⁶ These parametric conditions are achieved by using argon gas with a reference diameter $dg=4.17×10−10$ m and molecular mass $mg=6.63×10−26$ kg as the working fluid. The gas number densities of 2.45 × 10²⁵ and 2.45 × 10²⁴/m³ are used for the Kn = 0.35 and 0.035 cases, respectively, while the gas free-stream temperature, as well as the particulate surface temperature, is maintained at 300 K. Since the study entails free-stream gas flow conditions in subsonic as well as supersonic regimes, special subsonic boundaries are necessary for low-speed gas flow conditions in DSMC. Typical DSMC subsonic boundaries combine extrapolated flow properties from within the domain to determining the upstream characteristics required at the subsonic boundaries. However, this method incurs high scattering, which would result in numerically fluctuating outflow at the domain boundaries. We instead use far-field half-Maxwellian inflow boundary conditions, wherein particles fluxing through all the boundaries in the simulation domain are obtained from the Maxwell–Boltzmann distribution with inflow (free-stream) gas conditions. Extensive details regarding the implementation of the half-Maxwellian boundary condition are given in  Appendix. A three-dimensional cubical domain of side length “a” is used as the simulation domain in this study as shown in Fig. 2. The particulate geometry is placed at the geometric center of the domain as indicated by the representative fractal aggregate in the figure. The simulation domain edge lengths of a = 20 and 40 μm are used for the Kn = 0.035 and 0.35 cases, respectively, to reduce the effect of subsonic boundaries on the fluid-induced quantities on the particulate based on a separate sensitivity study performed by varying the edge length a of the simulation domain. Note that too large a value of “a” unnecessarily increases the computational load by increasing the total number of simulated DSMC gas particles in the domain.

The particle Reynolds number $Re=ρg,∞U∞dp/μ$ determined from the free-stream gas density $ρg,∞$⁠, free-stream gas velocity $U∞$⁠, monomer diameter d_p, and gas viscosity μ varies approximately from 5 to 280 for the Mach number–Knudsen number combinations used in this study. Since the length-scale d_p used in the definition of Kn and Re is the same, the $Kn−Re−Ma$ relation, $Kn=(Ma/Re)γπ/2$⁠, can be used to interrelate the Knudsen, Mach, and Reynolds numbers. The DSMC gas-phase evolves with a time step magnitude $Δt=5×10−11$ s, and an adaptive octree-based spatial discretization with five levels of refinement is used for all the simulations performed in this study. The refinement strategy generates approximately 100 000 and 180 000 DSMC collision cells in the slip and transitional regimes, respectively. A minimum of ten DSMC simulated particles are maintained within these cells throughout the simulation. The simulations are sampled over 20 000 time steps after the gas-particulate simulation system attains a steady-state in 20 000 simulation time steps to reduce the statistical noise associated with the particle-based DSMC method.

Before studying the mobility of complicated fractal geometries, the gas kinetic method is used to understand the mobility of an isolated spherical monomer of diameter d_p = 2 μm that will later be used to construct fractal aggregates. The fluid-induced drag force (along the flow direction) developed on the spherical particulate exposed to flows with Ma < 6.5 for the two specific Knudsen numbers Kn = 0.035 and 0.35 under consideration in this study is compared to the available semi-empirical correlations from the literature.

Figure 3(a) compares the coefficient of drag obtained from DSMC simulations for the intermediate Reynolds numbers (low Mach numbers) with Clift-Gauvin (CG)⁴⁰ and Loth's⁴¹ modified CG relation, which is given by the following two-part formulation:

The coefficients for the compressibility-dominated regime, Re > 45 is given as follows:

The coefficients in the rarefaction dominated regime, Re < 45 are given as follows:

where $s=Maγ/2$ represents the particle speed ratio. Note that Loth's model that accounts for the gas rarefaction and compressibility in the vicinity of the spherical particulate reduces to the CG model for $HM,CM,GM=1$⁠.

The variation of slip correction factor $CC=FStokes/FDSMC$ obtained from the DSMC calculations for different Knudsen number flows in the Stokes regime (⁠$Re≪1$⁠) is provided in Fig. 3(b). Here, $FStokes=3πμdpU∞$ and $FDSMC$ are the theoretical Stokes drag, and DSMC computed drag force, respectively. From the gas kinetic simulation perspective, argon gas viscosity μ is determined using the VHS formulation,³⁶ and free-stream gas velocities of $U∞=5$ and 1 m/s are used in high Kn > 0.5 and low Kn < 0.5 cases, respectively, to impose Stokes flow in the vicinity of the spherical particulate. The decrease in Stokes drag due to the reduced number of gas-particulate collisions in the rarefied or high Kn gas regimes captured by the DSMC simulations is in good agreement with the Loth's model in the low Re regimes that approach the semi-empirical formulations of Phillips et al.⁴² and Annis et al.⁴³ given as follows:

The high-fidelity DSMC results are in good agreement with the drag model proposed by Loth et al.⁴¹ as shown in Fig. 3 for the gas parametric space under consideration. The effect of gas rarefaction and compressibility on spherical particle drag predicted by Loth's model that modifies the generic CG drag correlation also is accurately captured by the DSMC method. Since there are no theoretical C_D formulations for these intermediate Reynolds numbers for variable Kn due to the nonlinearity effect of the gas-particulate momentum transfer,⁴¹ the two-part drag correlation by Loth et al.⁴¹ is used as the baseline case for the mobility parameters of the fractal aggregate to be discussed in Secs. IV–VII. Due to the statistical noise associated with the DSMC method at low Re or Stokes regime, a maximum difference of 15 % is observed from the theory. At higher Re, the difference is less than 5 %. Therefore, a similar accuracy can be expected for the mobility study of fractal aggregates in Sec. V B.

Fractal aggregates (FA) of complex geometry can be composed of monomers of varying diameters and geometries (spherical, spheroidal, ellipsoidal, etc.) glued together through various aggregation mechanisms. Simplified fractal-like aggregates are constructed from the arrangement of monodisperse spherical monomers of diameter $dp=2 μ$m in a three-dimensional space using a tunable particle cluster aggregation algorithm.⁴⁴ The FA thus constructed obeys the following scaling law:

where N is the number of monomers constituting the FA, the radius of gyration of the aggregate that describes the spatial distribution of the monomers is represented by R_gyr, the radius of the aggregating monomer is given by $rp=1 μ$m, and k_f and $Df<3$ represents the packing factor and the fractal dimension, respectively. A FA with dimension D_f closer to three has a spatial geometry that approaches a spherical geometry, while $Df∼1$ aggregates have an approximate chain-like construct. To generate fractal-like aggregates in computer simulations, researchers, such as Fillipov et al.⁴⁵ and Skorupski et al.,⁴⁶ have used stochastic aggregation methods. These methods can be classified into two categories: (1) particle–cluster aggregation (PC) where particles are attached to a growing cluster at every time step of aggregate generation and (2) cluster–cluster aggregation (CC) where smaller clusters aggregate to form larger ones every time step. In this work, PC aggregation is used to generate FA aggregates. Although PC aggregation has several disadvantages, such as directional bias, these methods are relatively efficient in generating small aggregates (⁠$N<100$⁠). Note that both these methods (PC and CC) are computationally very expensive and the computational load increases drastically with higher N and lower D_f.

The PC aggregation algorithm used in this study introduces a single spherical monomer or particle at every time step of the aggregation process. At each step of the aggregation process, the distance between the geometric center of the growing aggregate to the center of the new particle is specified by $|χ→|$ given in the following equation:

where N_p is the number of p-type particles or monomers already in the aggregate with the $p+1th$ particle being introduced at a radial location of $|χ→|$ with respect to the center of mass of the p-particle aggregate. The FA generated from the PC algorithm follows the scaling law given in Eq. (12). The center of mass $RCM→$ of the (p particle) aggregate assuming a uniform aggregate material density is given as follows:

where R_i is the center of mass of each individual monomer in the aggregation process. Figure 4 shows a representative step in the process of particle–cluster aggregation. The red particle labeled M in the figure is the $p+1th$ particle to be attached to the already existing p-particle aggregate. A random position on the sphere of radius χ computed from Eq. (13) is sampled in the Cartesian coordinate system using the method proposed by Marsaglia et al.⁴⁷ First, we select two variables X₁ and X₂ from a uniform distribution $[−1,1]$ such that $S=X12+X22<1$⁠. The random point R on the sphere is defined by the following formula:

where $î,ĵ$⁠, and $k̂$ are the mutually orthogonal basis vectors of the Cartesian coordinate system. To identify a successfully sticking or aggregating M particle position, a tolerance of ϵ is used such that the distance between the aggregating surface and the new position is less than $ε=0.1 μ$m.

The tunable particle-clustering aggregation code was used to generate a dense FA of $Df∼2.2$ and an open FA of $Df∼1.8$ with N = 10 primary particles or monomers each. The monomers forming both the aggregates were assumed to have spherical geometry and uniform diameters of 2 μm. It was observed that the generation of the open aggregate required substantially longer computer time (⁠$∼45$ min) as compared to denser aggregates. This observation can be attributed to the increased χ for the open aggregate, making the random sampling step in the clustering algorithm, Eq. (15), computationally expensive. The dense and open fractal aggregates used in the study, generated from the tunable PC aggregation code, are shown in Fig. 5.

The random arrangement of the primary particles in three-dimensional space creates interstitial spaces or voids in the aggregate, as is evident in Fig. 5, and it can be seen that the open aggregate is more porous than the dense aggregate. The flow characteristics of the gas in the vicinity of the aggregate will depend on the permeability of the FA and will be an important factor in determining the mobility parameter of the aggregates.

Since the FAs are directionally asymmetric, three mutually perpendicular orientations are used in the gas kinetic study to determine the general, directionally independent mobility parameters for the aggregates. The three orientations for the aggregates with respect to the flow direction are indicated in Fig. 5 by the thick arrows which also indicate the respective incoming flow direction. The projected area of the dense aggregate with respect to the flow direction for orientations 1, 2, and 3 are 20.9, 21.2, and 24.5 μm², respectively. Similarly, for the open FA along 1, 2, and 3 orientations, the projected area is 25.1, 29.06, and 15.7 μm², respectively. The dense aggregate has less variation between the projected areas at different orientations due to its dense packing of monomers. The branched or chained arrangement of primary particles brings about the larger variation of projected areas at different orientations for the open FA.

The complicated flow pattern and the resultant flow-induced surface force for a selected set of simulations of the Re and Kn number space for a fractal aggregate are shown in Fig. 6. The gas streamlines that decelerate through momentum transfer with the FA surfaces are classified as either an exterior streamline, which moves along the outer surface of the aggregate (indicated by e), or an interior streamline, which percolates through the inner voids of the aggregate (indicated by i). For the low Re flows shown in Figs. 6(a) and 6(c), the exterior and interior streamlines decelerate at the same rate due to the dominant viscous resistance at these low gas-speeds, while at higher Re, shown in Figs. 6(b) and 6(d), the effect of the fluid inertia is comparable to their viscous resistance. The tortuous gas streamlines moving through the interior voids of the FA decelerates to a much lower gas velocity as opposed to the exterior gas streamlines due to a larger number of gas–surface collisions that increase the momentum transfer between the gas and aggregate surface. In contrast, the gas molecules in the exterior streamlines constantly replenish their momentum through gas–gas collisions with the high momentum free-stream gas molecules. Since the momentum loss through gas–surface collisions is readily replenished by gas–gas collisions, these streamlines do not decelerate to the same extent as the interior ones.

The surface distribution of gas-induced drag force on the aggregate in the gas flow direction given by Eq. (2) is shown in Fig. 6 for gas Knudsen numbers of 0.035 and 0.35, respectively. For low Reynolds number or low speed flows in both the Knudsen numbers, the gas-induced force is distributed uniformly along the FA surface due to the slow viscous flow that transfers momentum uniformly along the particulate surface. At higher gas free-stream velocities or higher Re, the high-speed, high-momentum gas molecules directly impact the flow facing side of the aggregate, while the surfaces away from this side are exposed to the slower, low-momentum gas molecules. Due to this inertial effect of the gas, the surface force distribution shows a higher magnitude on the flow-facing side compared to the interior surfaces as well as surfaces in the lee of the flow.

In the compressible flow regimes, (Ma > 1), the gas molecules reflecting from the FA surface collide with the incoming molecules to develop a shock discontinuity. For the dense FA, due to the proximity of the monomers in the arrangement, the diffuse (bow) shock generated at the flow facing side of the aggregate is observed to be a collective, single entity as shown in Fig. 7(a), while for the open aggregate, individual shocks are observed on the branched ends of the aggregate as indicated in Fig. 7(b).

The mobility parameters of the FAs obtained by integrating the flow-induced surface forces and moments generated on the FAs using Eqs. (2) and (3) in the three-dimensional Cartesian coordinate system are resolved along the gas flow direction and its perpendicular by defining a Lagrangian coordinate system for the FAs. The coordinate system depicted in Fig. 2 resolves the forces and moments generated on the body into axial (P1) and lateral (P2) directions as follows:

where F_x, F_y, and F_z are the X, Y, and Z components of the flow induced force given by Eq. (2). M_x, M_y, and M_z are the X, Y, and Z components of the pitching torque given by Eq. (3), and the axial and lateral force components in the above equation are the drag force, F_d and lift force, F_l on the aggregate. These mobility parameters are further classified as a primary mobility parameter, e.g., drag force (F_d), and secondary mobility parameters, e.g., Moments (⁠$MAxial$ and $MLateral$⁠) and lateral force (F_l).⁴⁸ It is to be noted that the direction of the lift and lateral pitching moment given by Eq. (16) may not coincide. Additionally, to determine the relative contribution of viscous and pressure force on the FA drag and lift, the shear and pressure components of the surface force distribution are integrated separately along the P1 and P2 directions, respectively, while post-processing the DSMC simulation results.

The gas-induced drag force experienced by the two FAs at different orientations is shown in Figs. 8(a) and 8(b) for the slip and transitional regimes respectively through the normalized parameter $p=2Fd/ρ∞U∞2AmonoCD$⁠. Here, the drag force on the FA obtained from DSMC, F_d, is normalized with the drag on a single spherical monomer of cross-sectional area $Amono$ computed using the semi-empirical drag model by Loth et al.⁴¹ given in Eq. (5), and $ρ∞$ and $U∞$ are the free-stream gas density and flow velocity, respectively. Note that the normalization of F_d using the DSMC computed force on a spherical monomer given in Sec. III differs by < 4% and 8.4 % for the slip and transitional regimes from the analytic formula-based normalization, respectively. Based on the approximate logarithmic trend of the parameter p with Ma, a non-linear model is proposed for the normalized drag parameter p for the drag force averaged over orientations, O1, O2, and O3 as follows:

where the non-linear fit coefficients A₁ and B₁ are obtained using non-linear regression and are tabulated in Table I for the two aggregates in the slip and transitional gas regime. The table also includes the variance (⁠$σp2$⁠) of the DSMC computed data points of p parameter for the ten Mach numbers with respect to the model fit. Note that the variance of DSMC calculated p for a spherical monomer in the slip and transitional regimes (22.74 and 17.96, respectively) are much higher than the corresponding variances with respect to the model fit, suggesting that the proposed model accurately captures the fractal aggregate drag compared to a sphere.

In general, p increases with increasing Ma for all the cases, suggesting that the drag experienced by the fractal aggregates increases much more rapidly with Mach number compared to a sphere as indicated in Figs. 8(a) and 8(b) for slip and transitional regimes, respectively. This increase may be attributed to the additional viscous drag experienced by these aggregates at high Mach numbers that does not occur in the case of spheres in the same flow regime.⁴⁹ The non-linear fits for the drag parameter p are higher for the open aggregate compared to the dense aggregate in the high Ma regimes due to its increased average projected area. The rate of increase in p with Ma is approximately captured by the fit parameter $(A1·B1)$⁠. The parameters indicate that in the slip regime, both the open and dense aggregates show a steep increase in p at lower values of Ma while in the transitional regime, p increases gradually for the open aggregate. The orientation-averaged contribution of viscous drag shown in Fig. 8(c) at different Ma indicates that its contribution is higher at lower flow speeds (⁠$Ma→0$⁠). The attached flow on the aggregates [e.g., Fig. 6(a)] manifests as a high viscous contribution for particles in the slip regime. The high slip velocity at the aggregate surface due to the decrease in gas–surface collisions in the transitional regime decreases the percentage contribution of viscous drag in this regime for low Ma. In the low Ma regimes (Ma < 0.3) where the viscous effects are more prominent, the drag force experienced by either FA does not vary significantly with orientation as shown in Figs. 8(a) and 8(b). The increase in viscous contribution with Ma pronounced in the transitional regime (Kn = 0.35) can be attributed to the increase in viscous drag induced by the slow-moving gas through the interiors of the FAs behind the diffuse bow shock. For the higher Ma values, the viscous contribution is low for all the cases, and the surface-normal pressure drag is the major contributing factor. Since the pressure component of drag is a strong function of the projected area, the drag experienced by the FAs in this regime varies with orientation. For example, the drag value variance is approximately $52%$ for the dense FA at Ma = 6 compared to 12% at $Ma=0.2$⁠. In this regime, the open FA experiences a higher drag force than the dense FA due to the former's higher projected area compared to the latter.

At Kn = 0.035 and Ma = 6, the average drag force experienced by the open FA is 19.59% higher than the dense FA. Even though the drag force experienced by the FAs in the transitional gas inflow is lower than the slip regime due to less number of gas–surface interactions, the variation in the drag parameters at the lower and higher Mach numbers follows the same trend as that of the slip regime discussed above.

Next, the secondary mobility parameters, i.e., the lateral lift and pitching torque generated on the FAs, are presented. The variation of monomer-drag normalized lift parameter $q=2Fl/ρ∞U∞2AmonoCD$ is shown in Figs. 9(a) and 9(b) for the slip and transitional regimes, respectively. Here the lateral lift force experienced by the FA, F_l is normalized with the drag on a single spherical monomer of cross-sectional area $Amono$ computed again using the semi-empirical drag model by Loth et al.⁴¹ given in Eq. (5). Note that the normalization of F_l with the DSMC computed drag on a spherical monomer differs by < 4% and 8.4% for the slip and transitional regimes from the analytic formula-based normalization, respectively. Unlike the drag parameter p, the lift parameter q, which is an order of magnitude lower than the former, initially increases at lower Ma values, whereas at higher Mach numbers, it gradually decreases. Note that a similar order of magnitude difference in lateral force was observed by Chinnappan et al.²³ for prolate ellipsoidal particles exposed to a free-molecular gas flow. To capture this general trend, the following non-linear function is proposed for the orientation-averaged lift parameter as follows:

where the non-linear fit parameters A₂, B₂, and C₂ are determined from non-linear regression of the orientation-averaged DSMC lift force on the aggregate as given in Table I for the slip and transitional regimes. The table also shows the variance (⁠$σq2$⁠) of the DSMC computed lift parameter q with respect to the model fit. The numerator and denominator of the above function capture the initial increase and later decrease in the parameter q with Ma.

The decrease in lateral lift with respect to the monomer drag at higher Mach numbers may be attributed to the general decrease in pressure and viscous drag in the direction of lift (perpendicular to flow) as the flow perceives the aggregate as a spherical body at higher Ma. The rate of increase in q with Ma approximated as $(A2.B2)$ for low Ma indicates that q increases steeply in the slip regime for the open aggregate compared to the dense aggregate, while in the transitional regime, the rate is approximately the same. At high Ma numbers, the rate of decrease in q with Ma can be approximated as $(A2/C2)$ suggesting that q for the open aggregate decreases faster than for the dense aggregate. Figure 9(c) suggests that the viscous contribution to the lateral lift is significant at subsonic flow velocities and approximately 8 to 30 % higher than the viscous contribution for the particulate drag in the same regime. At higher Mach numbers, the viscous component is lower than at subsonic Mach numbers because the pressure component becomes the major contributor to the lateral lift. In this regime, the viscous contribution for the lateral lift F_l is approximately $10%$ lower than that of F_d, as shown in Figs. 9(c) and 8(c). Unlike Fig. 8(c) where the viscous component for F_d for Ma > 2 is higher for the transitional case, Fig. 9(c) suggests that the percentage viscous contribution to F_l is not significantly affected by the degree of gas rarefaction. Similar to the discussion about F_d, the effect of orientation on F_l is not significant at lower flow speeds, where the viscous drag is dominant as shown in Figs. 9(a) and 9(b) while at higher Ma, due to the dominant pressure contribution, the effect of orientation becomes evident. For example, the average variation of F_l for a $Ma∼0.3$ dense aggregate in the slip regime is approximately $11%$ while at Ma = 6, it is approximately $63%$⁠.

Figure 10 shows the variation of axial pitching torque $MAxial$ and lateral pitching torque $MLateral$ with Ma. We observe that the pitching torque generated along an axis in the axial direction, i.e., P1 is an order of magnitude lower than that of the pitching torque component generated in the lateral direction P2 regardless of Kn, orientation, or aggregate type by comparing Figs. 10(a) and 10(b) and 10(c) and 10(d). In general, the pitching torque acting on the FAs increases with increasing Ma or flow velocity. The axial pitching torque also shows no significant variation with change in orientation, Kn or the aggregate dimension D_f. As previously observed for the primary drag and the lateral lift, at low gas velocities, the variation of lateral pitching torque between the different orientations is negligible due to the dominant viscous contribution compared to the pressure contribution, but at higher Ma, a significant variance is observed between the lateral torque generated at different orientations. We also observe that the average lateral pitching torque generated on the open FA is approximately $15%$ higher than the dense aggregate at higher Ma. This observation is due to the branched, fluffy construct of the open FA compared to the dense FA that generates unbalanced forces at larger distances from the center of mass of the aggregate generating higher moments.

To understand the effect of internal porosity on the mobility parameters of the fractal aggregates, the gas-kinetic study similar to Sec.V is performed on an ellipsoid of regular geometry approximating the structure of the FA. For this study, the open FA in orientation 2 (Ref. Fig. 5) is chosen and compared with flow over the minimum volume ellipsoid approximating the FA constructed using the Löwner–John algorithm, as shown in Fig. 11. The constructed ellipsoid has principle radii $(a,b,c)=(2.86,3.04,5.75) μ$m and is positioned at the same orientation as the original open FA with respect to the flow direction. The ellipsoid whose volume is ∼five times larger than the original aggregate is exposed to a slip gas at various inflow Mach numbers in the range of 0–6.5 corresponding to the parametric space used in Sec. V.

The surface drag force distribution on the approximate ellipsoidal particulate is shown in Fig. 12 where it can be seen that the surface force on the flow facing side progressively increases with increasing Ma, as was also observed in the case of fractal aggregates (Fig. 6). The net forces experienced by the ellipsoid $F→elli$ is computed using Eq. (2). The drag and the lateral lift force experienced by the ellipsoidal particulate obtained by resolving the net force $F→elli$ in the flow-parallel and flow-perpendicular directions [Eq. (16)] is normalized with the corresponding forces generated on the FA as follows:

where $Fd,elli$ and $Fd,FA$ are the drag force experienced by the ellipsoidal particulate.

The variation of non-dimensional parameters $ηdrag$ and $ηlift$ with Mach number is shown in Fig. 13. In general, the parameters have magnitudes greater than one since the surface area of the ellipsoid is ∼1.45 times that of the corresponding FA, and it also has a larger projected area. The sensitivity of the drag and lift forces on the FA to the inner gas streamlines can be qualitatively inferred from the variation of these parameters. The approximate constant magnitude of $ηdrag$ across Ma indicates that the FA drag is less sensitive to the gas flow through the inner pores of the aggregate. In contrast, the lift parameter $ηlift$ has a steep increase at low Mach number, but at higher Mach numbers, it is more or less constant. The fractal lift force is sensitive to the inner streamlines in the low Mach number regimes, while at higher Mach numbers, their contribution is negligible and the flow perceives the aggregate as a single entity indicated by the plateauing $ηlift$ parameter at high Ma.

The mobility parameters obtained from the high-fidelity DSMC study are first used to develop a simple, approximate model for the mobility of small fractal aggregates that accounts for their complicated asymmetric shape. The following assumptions are used in developing the model. The dynamics of the FA and the gas are assumed to be one-way coupled, i.e., the gas phase is not affected by the presence of the aggregate. The forces on the FA are resolved in the axial and lateral directions F_d and F_l and computed using the drag (⁠$CD,FA$⁠) and lift (⁠$CL,FA$⁠), respectively. The two components of force F_d and F_l are assumed to be independent of each other or de-coupled as their force magnitudes obtained from the kinetic study are an order of magnitude different. Since the rotational moment or the pitching torque on the aggregate is small, we assume that the effect of particle rotation is negligible.

Based on the above assumptions, the orientation-averaged coefficient of drag $CD,FA$ and lift $CD,FA$ is formulated as follows:

where p and q are the dimensionless scaling factors given by Eqs. (17) and (18), respectively, and were obtained from the orientation-averaged non-linear least-square fits for the drag and lift for the aggregates considered in this study at different Kn, and C_D is the coefficient of drag on a single sphere given by Loth et al.⁴¹ given in Eq. (5). Again, p and q account for the orientation-averaged shape effect of the FA as a function of Ma.

Next, we use the approximate orientation-averaged model developed for the aggregate to study a two-dimensional Rankine gas vortex-particle system, which is common in multiphase flow research.^50–52 The flow profile of a Rankine vortex is given as follows:

where $Uθ$ is the tangential velocity of the vortex (U_r = 0). Λ = 10 m²/s is the flow circulation and r is a radial distance from the center of the vortex. The inner core of the vortex has a radius $R=0.25$ m where the tangential gas velocity increases linearly with radial distance r. The solid particulate is placed at a distance of R_p = 0.1 m from the center of the vortex with zero initial velocity. The entire particulate is assumed to have a homogeneous material density of 2500 kg/m³. The evolution of the particle phase is studied by integrating Newton's equation of motion in a one-way coupled fashion. The momentum coupling parameters in the drag and lift direction is given by Eq. 20. The force of drag $F→drag$ and lift $F→lift$ experienced by the fractal aggregate is recast as follows:

where $V→rel$ is the relative velocity of the particle with respect to the gas at the particle center-of-mass location. The lateral lift force is modeled as a stochastic term to account for the orientation effects of aggregates by choosing a random vector direction $ϕ$ with respect to the x-direction through a uniform random number generator RAND. $î$ and $ĵ$ are the unit vectors along the x and y directions. Note that the drag force term in this model is independent of the lateral lift or lateral motion of the FA and is non-stochastic. The Euler's forward integration (first-order accurate) is used to generate the particle trajectory and the evolution characteristics. The control parameter (CFL number = U_p$Δtp$/ $Δxc$⁠) for the particle phase is maintained $<0.10$ by choosing a particle time interval value (⁠$Δtp$⁠) such that the particle of velocity U_p does not move more than a single fluid cell $Δxc$ at any time step in the simulation.

The evolution of a volume equivalent spherical particle in the slip and transition regime is shown in Figs. 14(a) and 14(d), respectively. Note that the time axis is plotted as a linear axis perpendicular to the $Uθ$ contour for the Rankine vortex. The spherical particles in both regimes show negligible radial displacement in 3 s of system evolution. In contrast, the fractal aggregates spiral out of the vortex core due to the increased drag force it experiences using the approximate non-linear drag parameter p. The evolution of the dense fractal aggregate in the slip and transition regime shown in Figs. 14(b) and 14(e), respectively, indicates that the dense aggregate achieves a higher average radial velocity of 0.196 m/s in the transitional regime as compared to 0.145 m/s in the slip regime. Similarly, for the open aggregate, the particulate achieves an average velocity of 0.136 and 0.232 m/s in the slip and transitional regime, respectively. It can be seen that the parachuting effect of the dense and the open fractal aggregates due to their shape effects is captured by the approximate model. The influence of lateral lift was not found to be significant in the evolution of these fractal aggregate trajectories as the rate of increase in lateral lift parameter q is much less than the drag parameter p for these conditions.

DSMC method was used to study gas flows in the vicinity of small fractal aggregates to understand their mobility parameters, such as drag, lateral lift, and pitching moment. Although the mobility of such asymmetric aggregates has been studied in the low Reynolds number (Stokes regime) and free-molecular regimes, their mobility is not understood well in the intermediate Reynolds number–Knudsen number regimes, to the best of our knowledge. This study focused on a parametric space of flow Mach number $Ma∈(0.1,6)$ for two Knudsen numbers Kn = 0.035 and 0.35 corresponding to slip and transitional gas regimes. The DSMC method that captures the compressibility and rarefaction effects dominant in this parametric space was used to study the mobility of an open aggregate of dimension $Df=1.8$ and a dense particulate of $Df=2.2$⁠. To understand the variability of the mobility parameters with aggregate asymmetry, three different orientations of the FAs with respect to the flow direction were considered in the gas-kinetic study. The simulations accurately resolved the gas molecule–particulate surface interactions on a molecular level to precisely capture the gas-particle momentum transfer by the interior and exterior gas flows in the aggregate vicinity. The accurate prediction of these gas features is crucial in the determination of the orientation-averaged mobility of the particulates.

Several key results were obtained from the study. The complicated gas flow around the fractal aggregates generates drag, lift, and pitching torque on the aggregates due to its asymmetric construct. The orientation-averaged drag force on the two fractal aggregates was normalized using the drag force on a spherical monomer given by Loth's⁴¹ model to obtain a drag parameter p, which was found to be a non-linear monotonically increasing function of the gas free-stream Mach number for both the slip and the transitional regime. This indicates that the drag force on the aggregate increases at a much higher rate than that of a sphere with increasing Ma. It was also observed that the viscous contribution to the fractal aggregate drag increases with increasing Ma as opposed to a monotonically decreasing contribution in the case of spherical particles. This observation may be attributed to the additional momentum transfer between the gas and particulate through the slow-moving gas in the interior of the porous aggregates. A similar study of the non-dimensionalized lift parameter q suggests that q increases initially for low Ma < 2 but at higher gas compressibility or high Ma, the magnitude of q decreases for both the aggregates for both Kn numbers. The viscous contribution toward the lateral lift of the particulate was found to be insignificant at higher Mach numbers as opposed to their contribution to particulate drag. On the contrary, the pitching torque was found to be negligible with respect to the drag and lift timescale for these small aggregates since the monomers in these aggregates are closely packed to generate significant moments about the particulate center of mass. Comparison of the drag and lift forces generated on a FA with an approximate ellipsoid revealed that the lateral lift force on the FA is more sensitive to the inner gas streamlines as compared to drag in the low Mach number regimes.

The lift (q) and drag (p) parameters were used in a non-linear Mach number-based correction factor to the Loth's⁴¹ model for a spherical monomer to approximately model the mobility of these small fractal aggregates in the parametric space. The model was shown to be useful in large length-scale Eulerian–Lagrangian^53,54 simulations to study the transport of such particulates of complicated geometries in convoluted, complex flow regimes. The model was used to study the evolution of these particulates in a simple, canonical Rankine vortex in a one-way coupled fashion. The results illustrate the disparity between the trajectories of fractal aggregates and volume equivalent spherical particles. The spherical particles tended to relax readily with the high-gradient flow of the canonical vortex while the fractal aggregates spiraled out of the inner core of the vortex. Since most practical multiphase systems have irregular particulates formed by aggregation, the state-of-the-art large-scale Eulerian–Lagrangian simulations require such models to accurately predict the evolution of the system. Future work will target expanding the parametric space to understand the effect of highly rarefied gas regimes on particulate mobility and to develop a correction scheme generic across different levels of gas rarefaction.

This work was sponsored by the Department of the Defense, Defense Threat Reduction Agency under Award No. HDTRA1-20-2-0001. The content of the information does not necessarily reflect the position or the policy of the federal government, and no official endorsement should be inferred. This work also used the Comet-GPU supercomputing resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE) TACC.

The authors have no conflicts to disclose.

A.V.M: Conceptualization, methodology, software, validation, formal analysis, writing—original draft, and visualization. D.A.L: Methodology, writing—review and editing, supervision, project administration, and funding acquisition.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The number of particles introduced at the boundaries, η, is calculated as follows:

where n is the number density of the gas in the free-stream corresponding to 1 atm gas pressure, $=2.6×1025$$/$ m³. The thermal velocity corresponding to the free-stream gas temperature is represented by V_mp and the speed ratio by $s=Vbulk/Vmp$⁠. The angle between the surface representing the boundary and the stream direction is indicated by θ

The incoming particles have velocities u such that the thermal velocity, $uth=u−Vbulk$ is sampled from the space $[−Vbulk,3Vmp]$⁠. An acceptance rejection probability given by

is used to decide if the gas computational particle velocity is acceptable in the simulation domain. A different velocity is sampled out of the space in case the first velocity was rejected.

If the boundaries are located sufficiently far from the test aggregate location, where the streamlines have recovered to their free-stream values, these boundary conditions should physically represent the flow about the particle and numerically, the solution should be independent of the boundary locations. The same condition has been used by Zhang et al.³⁵ to study the scalar friction factor for fractal aggregates. The researchers were successful in studying drag coefficient for fractal aggregates (along the flow direction) for different Knudsen numbers (Kn < 0.3).