Numerical study on the influence of outlet diameters on the discharge characteristics of dual-outlet silos

The gas phase model consisted of governing equations and constitutive equations. The governing equations for the gas phase include the continuity equation and the momentum equation, which are presented as equations (1) and (2), respectively. (1) ∂ ∂ t ( α g ρ g ) + ∇ ⋅ ( α g ρ g v g ) = 0 \frac{\partial }{\partial t}({\alpha }_{\text{g}}{\rho }_{\text{g}})+\nabla \cdot ({\alpha }_{\text{g}}{\rho }_{\text{g}}{{\bf{v}}}_{\text{g}})=0 where α _g is the gas volume fraction, v _g is the gas velocity vector, and ρ _g is the gas density. (2) ∂ ∂ t ( α g ρ g v g ) + ∇ ⋅ ( α g ρ g v g v g ) = − α g ∇ P g + ∇ ⋅ τ ¯ ¯ g − β ( v g − v s ) + α g ρ g g \frac{\partial }{\partial t}({\alpha }_{\text{g}}{\rho }_{\text{g}}{{\bf{v}}}_{{\bf{g}}})+\nabla \cdot ({\alpha }_{\text{g}}{\rho }_{\text{g}}{{\bf{v}}}_{{\bf{g}}}{{\bf{v}}}_{{\bf{g}}})=-{\alpha }_{\text{g}}\nabla {P}_{\text{g}}+\nabla \cdot {\overline{\overline{\tau }}}_{\text{g}}-\beta ({{\bf{v}}}_{{\bf{g}}}-{{\bf{v}}}_{\text{s}})\hspace{11.5em}+{\alpha }_{\text{g}}{\rho }_{\text{g}}{\bf{g}} where P _g is the gas pressure, τ ¯ ¯ g {\overline{\overline{\tau }}}_{\text{g}} is the gas phase stress tensor, β is the interphase momentum exchange coefficient, ν _s is the solid velocity vector, and g is the gravitational acceleration. On the right-hand side of the equation, the terms represent the pressure gradient, the divergence of the viscous stress, the interphase interaction force (e.g., drag), and the gravitational source term.

Considering that the silo discharge process exhibited typical characteristics of slow dense flow, the gas phase was modeled as an incompressible Newtonian fluid, with turbulence effects being neglected. Accordingly, the constitutive equation for the gas phase stress tensor can be expressed as follows: (3) τ ¯ ¯ g = 2 α g μ g D ¯ ¯ g − 2 3 α g μ g t r ( D ¯ ¯ g ) I ¯ ¯ {\overline{\overline{\tau }}}_{\text{g}}=2{\alpha }_{\text{g}}{\mu }_{\text{g}}{\overline{\overline{D}}}_{\text{g}}-\frac{2}{3}{\alpha }_{\text{g}}{\mu }_{\text{g}}tr({\overline{\overline{D}}}_{\text{g}})\overline{\overline{I}} where μ _g denotes the gas dynamic viscosity and D ¯ ¯ g {\overline{\overline{D}}}_{\text{g}} is the gas-phase rate of strain tensor, which can be calculated by (4) D ¯ ¯ g = 1 2 [ ( ∇ v g ) + ( ∇ v g ) T ] {\overline{\overline{D}}}_{\text{g}}=\frac{1}{2}{[}(\nabla {{\bf{v}}}_{\text{g}})+{(\nabla {{\bf{v}}}_{\text{g}})}^{T}]

The governing equations for the solid phase (i.e., the particulate phase) included the continuity equation, the momentum equation, and the granular temperature transport equation [32], which are expressed as equations (5)–(7), respectively. (5) ∂ ∂ t ( α s ρ s ) + ∇ ⋅ ( α s ρ s v s ) = 0 \frac{\partial }{\partial t}({\alpha }_{\text{s}}{\rho }_{\text{s}})+\nabla \cdot ({\alpha }_{\text{s}}{\rho }_{\text{s}}{{\bf{v}}}_{\text{s}})=0 where α _s is the solid volume fraction and ρ _s is the solid density. (6) ∂ ∂ t ( α s ρ s v s ) + ∇ ⋅ ( α s ρ s v s v s ) = − α s ∇ P g − ∇ P s + ∇ ⋅ τ ¯ ¯ s + β ( v g − v s ) + α s ρ s g \frac{\partial }{\partial t}({\alpha }_{\text{s}}{\rho }_{\text{s}}{{\bf{v}}}_{\text{s}})+\nabla \cdot ({\alpha }_{\text{s}}{\rho }_{\text{s}}{{\bf{v}}}_{\text{s}}{{\bf{v}}}_{\text{s}})=-{\alpha }_{\text{s}}\nabla {P}_{\text{g}}-\nabla {P}_{\text{s}}+\nabla \cdot {\overline{\overline{\tau }}}_{\text{s}}\hspace{11.2em}+\beta ({{\bf{v}}}_{\text{g}}-{{\bf{v}}}_{\text{s}})+{\alpha }_{\text{s}}{\rho }_{\text{s}}{\bf{g}} where P _s is the solid pressure and τ ¯ ¯ s {\overline{\overline{\tau }}}_{\text{s}} is the solid stress tensor. (7) 3 2 ∂ ∂ t ( α s ρ s Θ ) + ∇ ⋅ ( α s ρ s v s Θ ) = ( − P s,kc I ¯ ¯ + τ ¯ ¯ s,kc ) : ∇ v s − ∇ ⋅ ( κ s ∇ Θ ) − J coll − J vis \frac{3}{2}\left[\frac{\partial }{\partial t}({\alpha }_{\text{s}}{\rho }_{\text{s}}\Theta )+\nabla \cdot ({\alpha }_{\text{s}}{\rho }_{\text{s}}{{\boldsymbol{v}}}_{\text{s}}\Theta )\right]=(-{P}_{\text{s,kc}}\overline{\overline{I}}+{\overline{\overline{\tau }}}_{\text{s,kc}}):\nabla {{\boldsymbol{v}}}_{\text{s}}-\nabla \cdot ({\kappa }_{\text{s}}\nabla \Theta )-{J}_{\text{coll}}-{J}_{\text{vis}}

On the right-hand side of equation (7), the four terms correspond respectively to the granular kinetic energy fluctuations induced by shear deformation, the spatial transport of granular temperature due to diffusion, the dissipation of kinetic energy due to inelastic collisions, and the viscous dissipation arising from interparticle friction and particle-fluid interactions. In this equation, Θ denotes the granular temperature, which characterizes the intensity of particle velocity fluctuations, P _s,kc represents the kinetic-collisional pressure, τ ¯ ¯ s,kc {\overline{\overline{\tau }}}_{\text{s,kc}} denotes the kinetic-collisional stress tensor, and κ _s represents the diffusion coefficient. The corresponding expressions for these terms are given as follows: (8) Θ = 1 3 v s 2 ¯ \Theta =\frac{1}{3}\overline{{{\boldsymbol{v}}}_{\text{s}}^{2}} (9) P s,kc = α s ρ s Θ + 2 ρ s ( 1 + e ) α s 2 g 0 Θ {P}_{\text{s,kc}}={\alpha }_{\text{s}}{\rho }_{\text{s}}\Theta +2{\rho }_{\text{s}}(1+e){\alpha }_{\text{s}}^{2}{g}_{0}\Theta (10) τ ¯ ¯ s,kc = 2 α s μ s,kc D ¯ ¯ s + α s λ s,kc − 2 3 μ s,kc t r ( D ¯ ¯ s ) I ¯ ¯ {\overline{\overline{\tau }}}_{\text{s,kc}}=2{\alpha }_{\text{s}}{\mu }_{\text{s,kc}}{\overline{\overline{D}}}_{\text{s}}+{\alpha }_{\text{s}}\left({\lambda }_{\text{s,kc}}-\frac{2}{3}{\mu }_{\text{s,kc}}\right)tr({\overline{\overline{D}}}_{\text{s}})\overline{\overline{I}} (11) κ s = 15 d s ρ s α s Θ π 4 ( 41 − 33 η ) 1 + 12 5 η 2 ( 4 η − 3 ) α s g 0 + 16 15 π ( 41 − 33 η ) η α s g 0 {\kappa }_{\text{s}}=\frac{15{d}_{\text{s}}{\rho }_{\text{s}}{\alpha }_{\text{s}}\sqrt{\Theta \pi }}{4(41-33\eta )}\left[1+\frac{12}{5}{\eta }^{2}(4\eta -3){\alpha }_{\text{s}}{g}_{0}+\frac{16}{15\pi }(41-33\eta )\eta {\alpha }_{\text{s}}{g}_{0}\right] where (12) η = 1 2 ( 1 + e ) \eta =\frac{1}{2}(1+e)

In equation (9), e is the restitution coefficient and g ₀ is the radial distribution function, which accounts for the increased probability of collisions between particles when the solid phase becomes dense. It is calculated as follows: (13) g 0 = 1 − α s α s, max 1 / 3 − 1 {g}_{0}={\left[1-{\left(\frac{{\alpha }_{\text{s}}}{{\alpha }_{\text{s,}\text{max}}}\right)}^{1/3}\right]}^{-1} where α _s,max represents the solid volume fraction at the maximum packing limit.

In equation (10), μ _s,kc denotes the shear viscosity of the solid phase, which consists of the kinetic and collisional contributions; λ _s,kc represents the bulk viscosity of the solid phase, characterizing the resistance of particles to compression and expansion; and D ¯ ¯ s {\overline{\overline{D}}}_{\text{s}} is the solid-phase rate of strain tensor. The detailed expressions of μ _s,kc, λ _s,kc, and D ¯ ¯ s {\overline{\overline{D}}}_{\text{s}} are given in equations (14) and (16), respectively. (14) μ s,kc = μ s,k + μ s,c {\mu }_{\text{s,kc}}={\mu }_{\text{s,k}}+{\mu }_{\text{s,c}} (15) λ s,kc = 4 3 α s 2 ρ s d s g 0 ( 1 + e ) Θ π {\lambda }_{\text{s,kc}}=\frac{4}{3}{\alpha }_{\text{s}}^{2}{\rho }_{\text{s}}{d}_{\text{s}}{g}_{0}(1+e)\sqrt{\frac{\Theta }{\pi }} (16) D ¯ ¯ s = 1 2 [ ( ∇ v s ) + ( ∇ v s ) T ] {\overline{\overline{D}}}_{\text{s}}=\frac{1}{2}{[}(\nabla {{\bf{v}}}_{\text{s}})+{(\nabla {{\bf{v}}}_{\text{s}})}^{T}]

In equation (14), μ _s,k and μ _s,c can be calculated from the following expressions: (17) μ s,k = 10 ρ s d s Θ π 96 α s ( 1 + e ) g 0 1 + 4 5 g 0 α s ( 1 + e ) 2 {\mu }_{\text{s,k}}=\frac{10{\rho }_{\text{s}}{d}_{\text{s}}\sqrt{\Theta \pi }}{96{\alpha }_{\text{s}}(1+e){g}_{0}}{\left[1+\frac{4}{5}{g}_{0}{\alpha }_{\text{s}}(1+e)\right]}^{2} (18) μ s,c = 4 5 α s ρ s d s g 0 ( 1 + e ) Θ π {\mu }_{\text{s,c}}=\frac{4}{5}{\alpha }_{\text{s}}{\rho }_{\text{s}}{d}_{\text{s}}{g}_{0}(1+e)\sqrt{\frac{\Theta }{\pi }}

In addition, particles in the silo remain in long-term contact and form a stable contact network under quasi-static or slow flow conditions, so that the bulk powder in the silo exhibits mechanical properties similar to those of a continuous solid. Therefore, the description of the solid-phase stress must account for the influence of interparticle friction. Accordingly, in this study, the solid-phase stress tensor τ ¯ ¯ s {\overline{\overline{\tau }}}_{\text{s}} and the solid pressure P _s were formulated as the superposition of kinetic and frictional contributions, as given by the following equations [33]: (19) τ ¯ ¯ s = τ ¯ ¯ s,kc + τ ¯ ¯ s,f {\overline{\overline{\tau }}}_{\text{s}}={\overline{\overline{\tau }}}_{\text{s,kc}}+{\overline{\overline{\tau }}}_{\text{s,f}} (20) P s = P s,kc + P s,f {P}_{\text{s}}={P}_{\text{s,kc}}+{P}_{\text{s,f}}

In equation (19), the frictional shear stress tensor τ ¯ ¯ s,f {\overline{\overline{\tau }}}_{\text{s,f}} is calculated using the following model [34,35]: (21) τ ¯ ¯ s,f = 2 α s μ s,f D ¯ ¯ s − 2 3 α s μ s,f t r ( D ¯ ¯ s ) I ¯ ¯ {\overline{\overline{\tau }}}_{\text{s,f}}=2{\alpha }_{\text{s}}{\mu }_{\text{s,f}}{\overline{\overline{D}}}_{\text{s}}-\frac{2}{3}{\alpha }_{\text{s}}{\mu }_{\text{s,f}}tr({\overline{\overline{D}}}_{\text{s}})\overline{\overline{I}} where μ _s,f denotes the frictional viscosity, which can be calculated from the following expression [35]: (22) μ s,f = 2 P s,f sin ϕ e 2 D ¯ ¯ s : D ¯ ¯ s + Θ / d s 2 {\mu }_{\text{s,f}}=\frac{\sqrt{2}{P}_{\text{s,f}}\hspace{.25em}\sin \hspace{.25em}{\phi }_{\text{e}}}{2\sqrt{{\overline{\overline{D}}}_{\text{s}}:{\overline{\overline{D}}}_{\text{s}}+\Theta /{d}_{\text{s}}^{2}}}

In equation (20), P _s,f denotes the frictional pressure, defined as the normal component of the frictional stress, and is calculated as follows: (23) P s,f = 0 α s < α s, min Fr ( α s − α s, min ) n ( α s, max − α s ) p α s ≥ α s, min {P}_{\text{s,f}}=\left\{\begin{array}{c}0\hspace{7.3em}{\alpha }_{\text{s}}\lt {\alpha }_{\text{s,}\text{min}}\\ \text{Fr}\frac{{({\alpha }_{\text{s}}-{\alpha }_{\text{s,}\text{min}})}^{n}}{{({\alpha }_{\text{s,}\text{max}}-{\alpha }_{\text{s}})}^{p}}\hspace{1em}{\alpha }_{\text{s}}\ge {\alpha }_{\text{s,}\text{min}}\text{}\end{array}\right. where Fr, n, and p are empirical coefficients related to material properties; α _s,min denotes the minimum solid volume fraction at which frictional effects become significant.

In this study, momentum exchange between particles and gas is primarily represented by the drag force. The drag force not only significantly influences the evolution of particle velocity but also acts as a coupling term in the momentum equations of both the solid and gas phases. The Gidaspow drag model [36] was adopted here to provide a unified description of drag in both dilute and dense regions. The Gidaspow model combines the features of the Wen and Yu’s model [37], suitable for dilute particle regions, and the Ergun’s model [38], applicable to dense particle regions. It employs a piecewise formulation based on the solid-phase volume fraction, as detailed below. (24) β = 3 4 C D α s α g ρ g ∣ v s − v g ∣ d s α g − 2.65 α s ≤ 0.2 150 α s 2 μ g ( 1 − α s ) d s 2 + 1 .75 ρ g α s ∣ v s − v g ∣ d s α s > 0.2 \beta =\left\{\begin{array}{c}\frac{3}{4}{\text{C}}_{D}\frac{{\alpha }_{\text{s}}{\alpha }_{\text{g}}{\rho }_{\text{g}}| {{\boldsymbol{v}}}_{\text{s}}-{{\boldsymbol{v}}}_{\text{g}}| }{{d}_{\text{s}}}{\alpha }_{\text{g}}^{-2.65}\hspace{4.7em}{\alpha }_{\text{s}}\le 0.2\\ 150\frac{{\alpha }_{\text{s}}^{2}{\mu }_{\text{g}}}{(1-{\alpha }_{\text{s}}){d}_{\text{s}}^{2}}+\text{1}\text{.75}\frac{{\rho }_{\text{g}}{\alpha }_{\text{s}}| {{\boldsymbol{v}}}_{\text{s}}-{{\boldsymbol{v}}}_{\text{g}}| }{{d}_{\text{s}}}\hspace{1em}{\alpha }_{\text{s}}\gt 0.2\text{}\end{array}\right. (25) C D = 24 α g Re p [ 1 + 0.15 ( α g Re p ) 0.687 ] {\text{C}}_{D}=\frac{24}{{\alpha }_{\text{g}}{\mathrm{Re}}_{\text{p}}}{[}1+0.15{({\alpha }_{\text{g}}{\mathrm{Re}}_{\text{p}})}^{0.687}] (26) Re p = d s ρ g ∣ v g − v s ∣ μ g {\mathrm{Re}}_{\text{p}}=\frac{{d}_{\text{s}}{\rho }_{\text{g}}| {{\boldsymbol{v}}}_{\text{g}}-{{\boldsymbol{v}}}_{\text{s}}| }{{\mu }_{\text{g}}}