Are Continuum Predictions of Clustering Chaotic?

TL;DR

This paper demonstrates that continuum models of gas-solid clustering in multiphase flows exhibit chaotic behavior, characterized by sensitive dependence on initial conditions and positive Lyapunov exponents, with chaos being dimensionally dependent.

Contribution

It provides evidence that continuum predictions of clustering are chaotic and links this chaos to secondary multidimensional instabilities rather than fundamental kinematic instability.

Findings

System exhibits sensitive dependence on initial conditions.

Largest Lyapunov exponent is positive, indicating chaos.

Chaotic behavior depends on the system's dimensionality.

Abstract

Gas-solid multiphase flows are prone to develop an instability known as clustering. Two-fluid models, which treat the particulate phase as a continuum, are known to reproduce the qualitative features of this instability, producing highly-dynamic, spatiotemporal patterns. However, it is unknown whether such simulations are truly aperiodic or a type of complex periodic behavior. By showing that the system possesses a sensitive dependence on initial conditions and a positive largest Lyapunov exponent, $\lambda_1 \approx 1/\tau$, we provide a tentative answer: continuum predictions of clustering are chaotic. We further demonstrate that the chaotic behavior is dimensionally dependent, a conclusion which unifies previous results and strongly suggests that the chaotic behavior is not a result of the fundamental kinematic instability, but of the secondary (inherently multidimensional) instability.

Equations

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot \varphi {\mathbf{U}}_s=0,$$

$$\frac{\partial (1-\varphi )}{\partial t}+\nabla \cdot (1-\varphi ){\mathbf{U}}_g=0,$$

$${\rho }_s\varphi \left(\frac{\partial {\mathbf{U}}_s}{\partial t}+{\mathbf{U}}_s\cdot \nabla {\mathbf{U}}_s\right)=-\varphi \nabla p_g-\nabla p_s$$

$$+\nabla \cdot {\sigma }_s-\beta \left({\mathbf{U}}_s-{\mathbf{U}}_g\right)+{\rho }_s\varphi \mathbf{g}$$

$${\rho }_g(1-\varphi )\left(\frac{\partial {\mathbf{U}}_g}{\partial t}+{\mathbf{U}}_g\cdot \nabla {\mathbf{U}}_g\right)=-(1-\varphi )\nabla p_g$$

$$+(1-\varphi )\nabla \cdot {\sigma }_g+\beta \left({\mathbf{U}}_s-{\mathbf{U}}_g\right)+{\rho }_g(1-\varphi )\mathbf{g}$$

$$\frac{3}{2}{\rho }_s\varphi \left(\frac{\partial T}{\partial t}+{\mathbf{U}}_s\cdot \nabla T\right)=\left({\sigma }_s-p_s{\mathbf{I}}_3\right):\nabla {\mathbf{U}}_s-\nabla \cdot \mathbf{q}$$

$$-\frac{3}{2}{\rho }_s\varphi T{\zeta }_1\left(\nabla \cdot {\mathbf{U}}_s\right)+\frac{3}{2}{\rho }_s\varphi \left(\xi -\frac{2\gamma }{m}T-{\zeta }_{0}T\right).$$

$${\lambda }_1^{\ast }=\underset{t^{\ast }\to \infty }{\lim }\underset{{\delta }_0\to 0}{\lim }\frac{1}{t^{\ast }}\ln \frac{\delta (t^{\ast })}{{\delta }_0},$$

$$\delta (t^{\ast })=\sqrt{\sum _{j=1}^{N_z}\sum _{k=1}^{N_y}\sum _{i=1}^{N_x}{\Vert {\tilde{\alpha }}_{i,j,k}(t^{\ast })-{\alpha }_{i,j,k}(t^{\ast })\Vert }^2},$$

$${\lambda }_1^{\ast }=\underset{n\to \infty }{\lim }\frac{1}{n\Delta t^{\ast }}\sum _{i=1}^n\ln \frac{{\delta }_n}{{\delta }_0}$$

$${\tilde{\alpha }}_{n+1}={\alpha }_n+({\delta }_0/{\delta }_n)({\tilde{\alpha }}_n-{\alpha }_n)$$

$$s_1$$

$$s_2$$

$$s_3$$

$$s_4$$

Full text

Granular matter is a collection of discrete, interacting solid particles which, like classic (molecular) matter, can be generally classified into one of three states Jaeger, Nagel, and Behringer (1996): i) under static conditions granular heaps or piles can sustain gravity-induced stress and behave like a solid Geng et al. (2001); ii) dense granular flows characterized by enduring and multi-particle contacts behave similarly to a fluid Forterre and Pouliquen (2008); and iii) rapid granular flows characterized by instantaneous contacts described as a granular gas Goldhirsch (2003). However, it is worth noting that all three granular states are only superficially similar to their molecular counterparts due to the dissipative nature of particle-particle contacts (inelasticity, friction, etc.) Jaeger, Nagel, and Behringer (1996); Goldhirsch (2003). In industrial operations, “rapid” (collision-dominated) granular flows are often encountered in devices in which the particles are fluidized by a gas Kunii and Levenspiel (1991); such gas-solid flows are the focus of this work. Rapid gas-solid flows are prone to an instability termed clustering in which particles tend to form spatially inhomogeneous patterns of high and low concentrations Fullmer and Hrenya (2017). Continuum or two-fluid models have long been known to be able to predict the qualitative nature of the clustering instability Agrawal et al. (2001) and more recent quantitative assessments have also shown promising results Mitrano et al. (2014); Fullmer and Hrenya (2016); Fullmer et al. (tted) However, it is yet unknown whether or not such predictions are chaotic. We take a first step in answering this question by simulating fluidization in an unbounded domain and calculating a positive largest Lyapunov exponent, thereby indicating that continuum predictions of clustering are in fact chaotic. Further, we show that chaotic behavior may be reduced to periodic behavior by constraining the dimensionality of the system.

I Introduction

Physicists and engineers have long been attracted to granular matter due to the rich variety of patterns and complex behavior they are prone to display Hopkins and Louge (1991); Goldhirsch and Zanetti (1993); Goldfarb, Glasser, and Shinbrot (2002); Conway, Shinbrot, and Glasser (2004); Aranson and Tsimring (2006); Vinningland et al. (2007); Cheng et al. (2008); Christov, Ottino, and Lueptow (2010); Seiden and Thomas (2011); Shinbrot (2015). Granular and multiphase gas-solid flows are also of central importance to energy, chemical, petrochemical, pharmaceutical, food processing and other key industrial sectors Kunii and Levenspiel (1991); Cocco, Karri, and Knowlton (2014). In this work, we are interested in rapid gas-solid flows, a regime characterized by nearly instantaneous and binary collisions, large Stokes numbers, $St\sim\rho_{s}Ud_{p}/\mu_{g}$, and relatively dilute solids concentrations, which commonly occur in the freeboards or risers of circulating fluidized beds, pneumatic conveying systems, etc.

Industrial particulate systems are often characterized by a large separation of scales van der Hoef et al. (2008); Li et al. (2013); Cocco et al. (2017) necessitating the use of continuum or two-fluid models (TFM) to make full, system-scale predictions in most cases. Unlike molecular dynamics (MD) or discrete element methods (DEM), which track the motion and collisions of every particle, continuum models only resolve average properties, e.g., solids volume concentration, mean velocity and granular temperature – a measure of the fluctuating kinetic energy of the solids phase. An analogy between particles and molecules is often employed using a kinetic theory (KT) approach to derive and constitute continuum models Gidaspow (1994); Brilliantov and Pöschel (2004); Pannala, Syamlal, and O’Brien (2010). We use the nomenclature KT-TFM here to indicate the multiphase nature of the continuum model, rather than the more common KTGF (kinetic theory of granular flows), which is often used interchangeably for granular flows (no interstitial phase) and multiphase flows.

The behavior of interest in rapid gas-solid flow is the clustering instability Fullmer and Hrenya (2017). Through the dissipation of granular energy and/or a mean relative motion between the phases, initially uniform distributions of particles tend to group together to form inhomogeneous distributions of dense clusters which are persistent yet dynamic in nature. With sufficient grid resolution, KT-TFM simulations are able to qualitatively predict the clustering behavior Agrawal et al. (2001). A growing body of work also suggests that continuum model predictions are also quantitatively accurate in the prediction of clusters Mitrano et al. (2014); Fullmer and Hrenya (2016); Fullmer et al. (tted).

It is reasonable to expect that direct numerical simulations (DNS) Yin et al. (2013) and CFD-DEM simulations Radl and Sundaresan (2014); Capecelatro, Desjardins, and Fox (2015), which resolve the particle-scale dynamics, should be chaotic as elastic Dellago and Posch (1997) and inelastic McNamara and Mareschal (2001) hard-sphere MD simulations are known be chaotic, even in the absence of clustering. For continuum predictions, however, the picture is not so clear. Several observations suggest that continuum predictions of clustering may be chaotic: i) visual inspection of the instantaneous patterns do not appear regular enough to be periodic Agrawal et al. (2001); Fullmer and Hrenya (2016); ii) continuum simulations of a wall bounded liquid-solid fluidized bed have been shown to exhibit a continuous solids concentration spectra Gevrin, Masbernat, and Simonin (2010); and iii) experimental measurements of gas-solid flows in risers are known to be chaotic Marzocchella et al. (1997); Ji et al. (2000) (although, clearly, physical gas-solid flows inherently discrete). Other observations, however, appear to indicate that continuum predictions of clustering are not chaotic: i) a previous study of clustering in a dilute riser flow (wall bounded fluidization) with a KT-TFM only revealed periodic behavior Benyahia, Syamlal, and O’Brien (2007); ii) it is possible that chaos in such simulations may only arise due to the on/off and/or very stiff (highly-nonlinear) behavior of the empirical models used to close the solid stresses near maximum packing Benyahia (2016); and iii) when caused by mean relative motion, the clustering instability in gas-solid flow is known to stem from a kinematic instability Fullmer and Hrenya (2017), a type of instability which is known to produce limit cycle behavior in gas-liquid fixed-flux TFMs Robinson et al. (2008); López de Bertodano et al. (2017).

Due to the competing suggestions noted above and in the absence of definitive evidence, it is unknown whether continuum predictions of clustering are chaotic or some type of complex limit cycle, e.g., a high-order n-torus. We seek to answer this outstanding question by studying KT-TFM predictions of clustering in gas-solid fluidization in a fully-periodic domain. We first study the nonlinear response of the system to small changes in the initial condition and then calculate the largest Lyapunov exponent from a pair of coupled simulations. Finally, we conclude by showing how the dimensionality of the system can affect the observed behavior, thereby resolving the aforementioned discrepancies.

II Two-Fluid Model

The KT-TFM used in this work is one of the most rigorous of such models to-date owing to the instantaneous particle-gas force included directly in the starting Enskog equation Garzó et al. (2012). For the sake of brevity, only the transport equations are included here:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The nine unknown variables to be solved for are the solids (volumetric) concentration, $\phi$, the solids- and gas-phase velocity vectors, ${\bf{U}}_{s}$ and ${\bf{U}}_{f}$, the gas pressure, $p_{g}$, and the granular temperature, $T$, a measure of the (isotropic) fluctuating kinetic energy in the solids phase. Solids-phase constitutive relations for the stress tensor, ${\bf{\sigma}}_{s}$, and granular heat flux, ${\bf{q}}$, and closures for the pressure, $p_{s}$, shear and bulk viscosity, granular conductivity, Dufour coefficient, and the first- and zeroth-order collisional cooling rates, $\zeta_{1}$, and $\zeta_{0}$, are derived from kinetic theory and can be found in the original work Garzó et al. (2012). Gas-solid interaction closures for the mean drag Beetstra, van der Hoef, and Kuipers (2007), $\beta$, the thermal drag, $\gamma$ and the neighbor effect $\xi$ are all derived from DNS data. (We note that the first-order thermal-Reynolds-number dependent term of the thermal drag model was re-fit to the original DNS data Wylie, Koch, and Ladd (2003) with a function that vanished in the zero concentration limit, $K(\phi)=\sqrt{0.3\phi}/(1-\phi)^{3.6}$.) The gas-phase shear viscosity, $\mu_{g}$, density, $\rho_{g}$ and the density of the constituent particles, $\rho_{s}$, are material properties and assumed constant, i.e., the gas-phase is treated as incompressible. The gravitational acceleration vector, ${\bf{g}}$ is also assumed constant and aligned in the vertical, $y$-dimension. Finally, it should be pointed out that although the instantaneous fluid-particle force is modeled with a Langevin equation Garzó et al. (2012), once the Enskog equation is averaged over the velocity distribution function the resulting KT-TFM is entirely deterministic.

The KT-TFM described above is solved numerically using the National Energy Technology Laboratory’s open-source MFiX code (https://mfix.netl.doe.gov/). MFiX is a finite-volume-based CFD code which solves the governing equations on a staggered grid with implicit Euler time advancement using a SIMPLE-type algorithm with variable time stepping. The Superbee flux-limiter is applied for all variable extrapolation.

The KT-TFM and the solution method described above, as well as the system and conditions described in Sec. III, are the same as those from a previous work where the primary focus was on validating the mean slip velocity against CFD-DEM data Fullmer and Hrenya (2016). However, several specific changes have been made to eliminate ”artificial” potential sources of chaos. No empirical solids (frictional) stress model is considered; even the pressure-like term used to limit the solids concentration below the maximum packing limit has been neglected here. Consequently, the radial distribution function (RDF) at contact of Carnahan and Starling Carnahan and Starling (1969) is used in place of the Ma and Ahmadi Ma and Ahmadi (1988) RDF. While there is little quantitative difference between the two RDFs for $\phi<0.5$, the latter does not have the appropriate asymptotic behavior as $\phi\to\phi_{\max}\approx 0.64$, only diverging in the unphysical limit of $\phi\to 1$. In the previous study Fullmer and Hrenya (2016) a dynamic flux renormalization, $j_{y}=\phi v_{s}+(1-\phi)v_{g}=0$, was applied iteratively to ensure that the frame of reference did not drift due to numerical precision, i.e., round-off errors . Except where noted, the renormalization has been removed here so that $j_{y}$ can, and does, drift from zero. Finally, the simulations reported here have been run in serial. Although serial computation increased the computational demands considerably (previous simulations Fullmer and Hrenya (2016) were parallelized over 36 cores), serial simulations remove suspicion that chaos may be attributed to domain decomposition errors.

III System and Conditions

Four nondimensional parameters determine the condition of the system: the Archimedes number, $Ar=\rho_{g}(\rho_{s}-\rho_{g})\left|{\bf{g}}\right|d_{p}^{3}/\mu_{g}^{2}$ where $d_{p}$ is the particle diameter, the density ratio $\rho^{*}=\rho_{s}/\rho_{g}$, the mean solids concentration $\left\langle\phi\right\rangle$ where the angle brackets indicate a volume average, and the restitution coefficient, $e$. While not appearing directly in Eqs. (1) - (5), the restitution coefficient describes the dissipation of energy due to particle collisions and appears in all of the solids-phase closures. In this KT-TFM Garzó et al. (2012), $e$ is treated as a constant material property and takes a value of $e=0.9$ here. The dimensional properties originally Radl and Sundaresan (2014) selected to represent typical fluid catalytic cracking particles in air correspond to $Ar=24.886$ and $\rho^{*}=1153.8$. Of the six mean concentrations considered previously Fullmer and Hrenya (2016), only the lowest, $\left\langle\phi\right\rangle=0.02$, is considered here to avoid unphysical predictions near maximum packing.

The gas-phase pressure is decomposed into a local fluctuating component and a constant linear component that balances the force of gravity. The domain size is $L_{x}^{*}=L_{z}^{*}=106.\bar{6}$ in the transverse directions and four times as tall in the vertical (streamwise) direction, $L_{y}^{*}=4L_{x}^{*}$. The simulations specify a relatively fine, uniform grid of $N_{x}\times N_{y}\times N_{z}=4\cdot 32^{3}$. In the previous study, simulations were carried out from a randomly perturbed initial state to a final condition at $t^{*}=76.8$. The final state is used as an initial condition in this work and the time is considered start at zero from this point in time. Here, space has been nondimensionalized by the particle diameter, $L^{*}=L/d_{p}$, and time is nondimensionalized by the viscous relaxation time, $t^{*}=t/\tau$, where $\tau=\rho_{s}d_{p}^{2}/18\mu_{g}$ is the viscous relaxation time.

IV Sensitive Dependence on Initial Conditions

As a preliminary assessment of whether or not the system might be chaotic, we begin by studying the response of the system to slight changes in the initial condition. Three simulations were considered. The reference case is unperturbed, i.e., uses the final condition from a previous study. Two perturbed cases were considered with nearly the same initial condition but with the granular temperature uniformly perturbed by $T_{0}\to T_{0}\pm\epsilon$ where $\epsilon=3.0\cdot 10^{-5}$ is a constant, the magnitude of which if further justified in Sec. V, and $T_{0}$ is the unperturbed initial granular temperature field. We choose to perturb only the granular temperature in an effort to disturb the system as innocuously as possible.

The evolution of solids concentration from the three slightly different initial states is displayed in Fig. 1 starting at $t^{*}$ = 7.68 and separated by the same amount. At the first two times (left half of Fig. 1), the three simulations are virtually indistinguishable. By $t^{*}=23.04$, the concentration in the top half of the domain of the perturbed case ${\tilde{T}_{0}}=T_{0}+\epsilon$ (top row) is clearly different than the reference case. At this time, the lower half of the domain remains relatively similar. The second perturbed case ${\tilde{T}_{0}}=T_{0}-\epsilon$ (bottom row) only displays minor deviations from the reference case. However by the time $t^{*}=30.72$, the three simulations are completely different from one another, resembling three entirely different realizations of the same flow. This simple demonstration shows that the KT-TFM (continuum) prediction of clustering in unbounded fluidization – in the absence of switches, functions containing singularities or highly-nonlinear empirical stress models – possesses a sensitive dependence on initial conditions, a strong hint that the flow may be chaotic.

V Largest Lyapunov Exponent

While no formal proof of chaos exists Sprott (2003), the Lyapunov spectrum remains one of the most generally accepted metrics in identifying chaotic behavior. Due to the size of this system, $N=32\times 128\times 32\times 9$ coupled ODEs, and its associated computational cost, we are only interested here in calculating the largest Lyapunov exponent (LLE), $\lambda_{1}^{*}$. A positive valued LLE indicates a chaotic system and, if positive, its magnitude determines the rate at which predictability of the system is lost. The LLE is defined by Cencini, Cecconi, and Vulpiani (2010)

[TABLE]

where $\delta_{0}$ is the initial perturbation and $\delta(t)$ is the evolution of the perturbation with time. By relaxing the inner limit of Eq. (6), which is meant to ensure linearity, one can obtain a rough estimate of the LLE using $\delta(t)$ from the two perturbed cases in Fig. 1. For this $N$-dimensional system, $\delta(t)$ is the L2-norm of the difference between the perturbed and reference trajectory summed over all primary variables, i.e.,

[TABLE]

where the tilde denotes the perturbed trajectory and $\alpha$ is the vector of unknown (solution) variables, ${\bf\alpha}=(\phi,u_{g},v_{g},w_{g},p_{g},u_{s},v_{s},w_{s},T)$.

The separation of the two initially perturbed cases relative to the reference case in Fig. 1 are compared quantitatively in Fig. 2. After an initial period of contraction ($t^{*}<1$), both cases diverge roughly exponentially from the reference case. At $t^{*}\approx 25$ and 30, both separations saturate when the three solutions have become essentially independent of one another, as shown in the fourth column of Fig. 1. It is unknown what causes the abrupt changes, notably just before $t^{*}=10$. It was verified that this was not an output or post-processing error nor did either of the calculations need to reduce the time step for convergence around this time. The divergence can be fit to a rate of approximately $0.2e^{0.64t^{*}}$ shown as a dashed line in Fig. 2. Although this analysis is admittedly imprecise since the separation becomes too large to be considered in a linear regime, it nonetheless indicates that the LLE may be in the vicinity of $\lambda_{1}^{*}\sim 0.64$.

A more precise calculation of the LLE requires maintaining the inner linear-constraining limit in Eq. (6). In practice, the LLE is calculated by iteratively re-normalizing the perturbed trajectory to the reference or fiducial trajectory iteratively so that the two remain close Wolf et al. (1985); Sprott (2003); Cencini, Cecconi, and Vulpiani (2010). Therefore, Eq. (6) becomes,

[TABLE]

where $\delta_{n}$ is the expansion between times $t_{n-1}^{*}$ and $t_{n}^{*}=t_{n-1}^{*}+\Delta t^{*}$. After every $\Delta t^{*}$ step, the perturbed trajectory is re-normalized from the reference trajectory,

[TABLE]

preserving the magnitude of the initial separation, $\delta_{0}$ and the orientation of the previous step. Therefore the separation from $t_{n-1}^{*}$ to $t_{n}^{*}$ always begins with the same magnitude of initial separation, $\delta_{0}$. To carry out these calculations, two separate simulations were coupled together, periodically exporting solution data from the reference to the perturbed simulation. The time step used for LLE calculation $\Delta t^{*}$ is not the actual code iterative time step, $dt^{*}$. As discussed in Sec. II, MFiX uses an adaptive time stepping to overcome non-convergence issues. It was determined that forcing the code time step to be small enough as to avoid non-convergence issues was prohibitively expensive for this calculation. Therefore, the code was allowed to freely adjust the time advancement step and was forced to sync to a constant $\Delta t^{*}=3.84\cdot 10^{-3}$. Typically, each $\Delta t^{*}$ step contains between two to five iterative time advancement steps. As in Sec. IV, the reference case is restarted from the final state of previous calculations Fullmer and Hrenya (2016) and the perturbed trajectory is initially offset by ${\tilde{T}_{0}}=T_{0}-\epsilon$ with $\epsilon=3.0\cdot 10^{-5}$ for an initial separation of $\delta_{0}=\sqrt{N}\epsilon/3$. The value of $\epsilon$ is chosen so that when normalized over all $N$ ODEs, the separation is two orders of magnitude larger than the residual convergence criteria of the iterative-implicit Euler time advancement scheme, here set at 10*-7*.

It takes approximately $t^{*}=25$ for the initial perturbation to propagate throughout the system and align itself in the direction of maximum expansion, marked by a rather obvious change in the behavior of $\delta_{n}$. Perhaps not coincidentally, this time corresponds to the saturation of the initial perturbation in Fig. 2. We begin collecting data for the calculation of the LLE at $t^{*}=38.4$ through the end of the simulation at $t^{*}=76.8$. While this coupled-simulation pair took over 67 days to complete, this yields only $n=10^{4}$ separations to average over, orders of magnitude smaller than other LLE measurements of spatio-temporal chaos Brummitt and Sprott (2009); Fullmer et al. (2014) and not sufficient to report a converged value with significant precision. However, in this study we are less concerned with the exact value of $\lambda_{1}$ than its sign and, to a lesser extent, its approximate value. Therefore, we choose to split the data into ten equally spaced intervals. Of the ten samples, the minimum, average, and maximum values are calculated to be $\lambda_{1}^{*}$ = 0.6981, 0.9767, and 1.1714, respectively. Again, the most important result is that we have demonstrated with reasonable certainty that $\lambda_{1}^{*}$ is positive, indicating a mean divergence of the perturbed trajectory from the fiducial trajectory and signifying chaotic behavior. The rough (nonlinear) approximation from Fig. 2 agrees relatively well with the $\lambda_{1}^{*}$ calculations, particularly the lower estimate. While approximate, the order of magnitude of $\lambda_{1}^{*}$ shows that the rate of predictability is lost inversely proportionally to the viscous relaxation time. It would be interesting to determine whether this finding holds in general, but is outside the scope of this exploratory study.

VI Discussion

The simple, albeit computationally intensive, measurement of a positive LLE indicates that KT-TFM (continuum) predictions of clustering are chaotic in nature. While this seems to agree intuitively with the physical picture of the clustering instability, it begs the question of what is different about this case than previous results that reported only limit cycle behavior. The answer appears to lie in the dimensionality of the system studied here. The previously mentioned cases for particulate gas-solid Benyahia, Syamlal, and O’Brien (2007) and bubbly gas-liquid Robinson et al. (2008); López de Bertodano et al. (2017) flows all considered a one-dimensional system. Here we provide a simple demonstration of the effect of dimensionality on clustering predictions by eliminating the $z$-dimensional dependence entirely and reducing the width of the system to $N_{x}$ = 1, 2 4 and 8 nodes. Note that the $N_{x}=1$ case is purely one-dimensional. With width of each node is the same as in previous cases. These four simulations were initialized at rest with a random, normal perturbation. The flux renormalization procedure from the previous study Fullmer and Hrenya (2016) was included so that limit cycle behavior could be identified without drift in the vertical flux. The reduced dimensionality simulations were run 10$\times$ longer than previously, up to times $t^{*}=768$.

For $N_{x}=1$, the system settles into a limit cycle pattern with three (relatively) uniformly spaced waves in the domain. Fig. 3 shows a snapshot of the concentration pattern (far left) and the nearly identical waves passing a probe point at ($x_{0}$,$y_{0}$) = ($d^{*}_{x}/2$, $L^{*}_{y}/2$), which appear incredibly qualitatively similar to the fixed-flux TFM bubbly void waves (e.g., see Fig. 4 of Ref. Robinson et al. (2008)). In order to visualize the limit cycle, consider the following four state variables:

[TABLE]

with state variables evaluated at ($x_{0}$,$y_{0}$). The ${\mathcal{H}}$ superscripts indicate the homogeneous solution Fullmer and Hrenya (2017) and $Re_{T}$ and $Re_{m}$ are the thermal- and mean-flow Reynolds numbers Garzó et al. (2012); Fullmer and Hrenya (2017). The linearly stable solution coincides with $s_{1}=s_{2}=s_{3}=s_{4}=0$, hence these four state variables provide a local measure of inhomogeneity. The $s_{1}-s_{2}-s_{3}$ phase-space colored by $s_{4}$ is shown in Fig. 4; the black sphere denotes the origin. For the case $N_{x}=1$, the data at ($x_{0}$,$y_{0}$) collapse on a limit cycle with no transverse features.

With just two transverse nodes, the one-dimensional limit cycle behavior is destroyed as illustrated in Fig. 4, although this case would need further study to say whether or not this case is chaotic or toroidal. However, the instantaneous wave pattern for $N_{x}=2$ displayed in Fig. 3 already hints at the distinction between periodic and chaotic behavior: i) the waves are not uniformly spaced, ii) the waves are not of equal amplitude and, most importantly iii) the lower of the three waves is not uniform in the transverse dimension. This behavior is more easily observed in Fig. 3 with $N_{x}$ = 4 and 8 nodes. In the (relatively) wider systems, transverse with spanning plugs also form, at least temporarily, analogous to the 1-D predictions. However, a secondary instability Batchelor (1993); Anderson, Sundaresan, and Jackson (1995); Glasser, Kevrekidis, and Sundaresan (1997); Glasser, Sundaresan, and Kevrekidis (1998); Duru and Guazzelli (2002); Sundaresan (2003) is then triggered which causes the plug to have a multi-dimensional funnel shaped structure which can be seen in Fig. 3. Ultimately, the plug spills through the funnel or, equivalently (depending on the readers preferred frame reference), the downstream bubble pushes through and breaks up the plug. The process can be qualitatively perceived as a translating Rayleigh-Taylor instability. In the phase-space of Fig. 4, the periodic behavior is destroyed entirely, the origin becomes obscured by the entangled dynamics and their magnitude grows away from the origin.

Although the root cause of clustering is a kinematic instability, that alone appears to be insufficient to cause chaotic behavior. Therefore, when the secondary instability is eliminated or suppressed through dimensional or domain-size reductions, only limit cycle behavior may be expected – at least for continuum predictions of multi-phase flow instabilities. It should further be noted that, while necessary, a multidimensional domain is not a sufficient condition for chaos. The governing conditions of the system studied here were sufficiently unstable for the multidimensional domain to be chaotic; other conditions, notably at higher $Ar$, may only yield periodic behavior even in multidimensional domains Gevrin, Masbernat, and Simonin (2010); Fullmer et al. (tted).

It is worth noting that the findings in this work carry some practical consequences for KT-TFM users. For one, chaotic behavior means that one cannot simply look at a single instance in time for quantitative information. Simulations need to be either averaged over a sufficiently long period of time or run as several realizations and ensemble averaged. This consequence is not very surprising, because, whether the dynamics were periodic or chaotic, it is universally known that most solutions are temporal and sufficiently complex to require averaging of measures and statistics. The second practical consequence is that slight differences between simulations – no matter how small, so long as the difference is representable by the precision of the floating point operations of the numerical solution – will eventually diverge into different instantaneous solutions. Differences may be something rather obvious like different initial conditions, different grids, different time steps, different numerical schemes, different tolerance limits, etc. However, solution differences may also result from unexpected sources such as compiling the same code on different machines, using different compilers or potentially something as simple as different compiler options. In closing, we remark that since bubbles in a dense emulsion seem to stem from the same instability mechanisms as clusters in a dilute suspension Glasser, Sundaresan, and Kevrekidis (1998), continuum predictions of bubbling gas-solid fluidized beds are also likely chaotic – although such a model would need to be augmented with several of the closures specifically neglected in this study.

Acknowledgements.

The authors are foremost grateful to Dr. Sofiane Benyahia; this study would have never been undertaken without his honest and forthright discussions with the MFiX community. We are grateful for the funding support provided by the U.S. Department of Energy under Grant No. DE-FE0026298.