Settling of inertial particles in turbulent Rayleigh-Bénard convection

The settling behavior of small inertial particles in turbulent convection is a fundamental problem across several disciplines, from geophysics to metallurgy. In a geophysical context, the settling of dense crystals controls the mode of solidification of magma chambers and planetary scale magma oceans, while rising of light bubbles of volatiles drives volcanic outgassing and the formation of primordial atmospheres. Motivated by these geophysical systems, we perform a systematic numerical study on the settling rate of particles in a rectangular two-dimensional Rayleigh-Bénard system with Rayleigh number up to ${10}^{12}$ and Prandtl number from 10 to 50. Under the idealized condition of spherically shaped particles with small Reynolds number, two limiting behaviors exist for the settling velocity. On the one hand, Stokes law applies to particles with small but finite response time, leading to a constant settling rate. On the other hand, particles with a vanishing response time are expected to settle at an exponential rate. Based on our simulations, we present a physical model that bridges the gap between the above limiting behaviors by describing the sedimentation of inertial particles as a random process with two key components: (i) the transport of particles from vigorously convecting regions into sluggish, low-velocity “piles” that naturally develop at the horizontal boundaries of the system, and (ii) the probability that particles escape such low-velocity regions without settling at their base. In addition, we identify four distinct settling regimes and analyze the horizontal distribution of sedimented particles. For two of these regimes settling is particularly slow and the distribution is strongly nonuniform, with dense particles being deposited preferentially below major clusters of upwellings. Finally, we apply our results to the crystallization of a magma ocean. Our prediction of the characteristic residence times is consistent with fractional crystallization, i.e., with the efficient separation of dense crystals from the residual lighter fluid. In absence of an efficient mechanism to reentrain settled particles, equilibrium crystallization appears possible only for particles with extremely small density contrasts.

Figure 1

Model parameter space describing particle dynamics. Red dashed and solid lines show three isolines of the terminal velocity $v_{\mathrm{t}}$. Colored symbols show three different simulation sets: A, B, and C (triangles, circles, and squares). The green line marks a constant value of required CPU time. We performed numerical simulations corresponding only to the sets B and C (see text for details). Pr is equal to 50.

Figure 2

(a) Volumetric averages of the root-mean-square velocities in simulation sets ${\mathrm{C}}_{50}^8, {\mathrm{C}}_{50}^{10}, {\mathrm{C}}_{50}^{12}, {\mathrm{C}}_{10}^8, {\mathrm{C}}_{10}^9$, and ${\mathrm{C}}_{10}^{10}$. The number in parentheses indicates the time-averaged value $u_{\mathrm{rms}}$ that is used later for computing the $v_{\mathrm{t}}/u_{\mathrm{rms}}$ ratio. (b) The corresponding Reynolds number, $U_{\mathrm{rms}}H/\nu$.

Figure 3

Settling curves obtained from the simulation set ${\mathrm{C}}_{50}^{10}$. All particle types are separated into four groups based on the $|v_{\mathrm{t}}|/u_{\mathrm{rms}}$ ratio, which is shown in color for each subplot. The $x$ axis represents the terminal distance $t{v}_{\mathrm{t}}$. Black and red solid lines represent the analytic solutions (14) and (15), respectively.

Figure 4

Selected settling curves from the simulation set ${\mathrm{xC}}_{50}^{10}$ that consists of particles with a high $|v_{\mathrm{t}}|/u_{\mathrm{rms}}$ ratio (2 to 50 for the depicted selection). Solid lines show the theoretical prediction (14), while colored points are obtained from the numerical simulation.

Figure 5

Temperature field in the simulation ${\mathrm{C}}_{50}^{10}$ at the time $t=12.9$. The temperature range is clipped for a better visibility of the up- and downwellings. Dots show particles with $0.3<v_{\mathrm{t}}/u_{\mathrm{rms}}<2.0$ that have not settled by the respective time. See Supplemental Material, Video S1 [61].

Figure 6

Same as Fig. 5, but at time $t=129.5$ and the depicted particles have smaller terminal velocities, $|v_{\mathrm{t}}|/u_{\mathrm{rms}}<0.3$. In green we show regions where thermal convection is slow, with velocities below $0.5u_{\mathrm{rms}}$. The thick green line separates the low-velocity piles from centers of large-scale convection rolls, which also show small velocities. See Supplemental Material, Video S2 [61].

Figure 7

Percentage of settled particles (green dots) compared with the two analytic predictions (14) and (15) (black and red lines, respectively). Blue stars show the percentage of particles that have entered low-velocity piles at least once by the respective time. The figure corresponds to time $t=50$.

Figure 8

(a) Settling probability $1-{\mathcal{P}}_{\mathrm{e}}$, where ${\mathcal{P}}_{\mathrm{e}}$ is the probability of escaping from a region with $\left|\mathit{u}\right|/u_{\mathrm{rms}}<p_{\mathrm{f}}$, with $p_{\mathrm{f}}=0.5$. Black symbols show $v_{\mathrm{t}}{t}_{\mathrm{c}}$ (or $-v_{\mathrm{t}}{t}_{\mathrm{c}}$ for light particles), where $t_{\mathrm{c}}$ is the average time between repeated entries into the low-velocity piles. (b) Ratio of $v_{\mathrm{t}}{t}_{\mathrm{c}}/(1-{\mathcal{P}}_{\mathrm{e}})$ for different values of the pile factor $p_{\mathrm{f}}$. The black dashed line shows the reference value 1 because for $v_{\mathrm{t}}{t}_{\mathrm{c}}=(1-{\mathcal{P}}_{\mathrm{e}})$ the solution (15) is obtained. The solid red line shows the skew normal distribution $f$ that is used in our analytic model for particle settling.

Figure 9

Solid lines show the temporal evolution of the average time between captures, $t_{\mathrm{c}}$, for selected $v_{\mathrm{t}}/u_{\mathrm{rms}}$ ratios. For a comparison, the characteristic time of a flow overturn is 5. The black line is equal to $max(1/0.3-t/2,0.5/0.3)$. Dashed lines show the corresponding percentages of settled particles.

Figure 10

Normalized residence time $F$ as a function of buoyancy ratio $\mathrm{\Lambda }$ and particle response time St. The residence time is computed with our analytic model (22), while squares, circles, and pentagons are colored according to the actual residence times obtained from the numerical simulations ${\mathrm{C}}_{50}^{10}, {\mathrm{xC}}_{50}^{10}$, and ${\mathrm{B}}_{50}^{10}$, respectively. The region of relatively slow settling (“slow belt” in the figure) occupies the area with $0.02\lesssim |v_{\mathrm{t}}|/u_{\mathrm{rms}}\lesssim 2.0$. Black triangles correspond to the parameters from Table 1. The color scale is clipped in order not to be dominated by the top right region, where settling occurs as if there was no convection.

Figure 11

Horizontal distribution of settled particles in the simulation set ${\mathrm{C}}_{50}^{10}$. The $v_{\mathrm{t}}/u_{\mathrm{rms}}$ ratio is depicted in color, and only the bilinear range, $0.3<|v_{\mathrm{t}}|/u_{\mathrm{rms}}<2.0$, is selected. (a) Vertically averaged vertical component of the fluid velocity, $u_z$. The quantity is further averaged over time, with the time window of integration covering complete sedimentation of the respective particles. (b) Heavy particles. (c) Light particles. Gray symbols indicate the flow direction of two major convection cells. The settled particles are sorted into 20 equally wide bins and the plotted value is normalized by the initial number of particles, i.e., divided by $N_0/20$.

Figure 12

As Fig. 11, only here we show the horizontal distribution of settled particles in the simulation set ${\mathrm{xC}}_{50}^{10}$.

Figure 13

As Fig. 12, only here we analyze the transitional regime, $0.02<|v_{\mathrm{t}}|/u_{\mathrm{rms}}<0.03$.

Figure 14

Same as Fig. 3, only here for the simulation sets ${\mathrm{B}}_{50}^{12}$ (left) and ${\mathrm{C}}_{50}^{12}$ (right). While ${\mathrm{B}}_{50}^{12}$ contains particle types with the modified density ratio $\beta$ very close to 1 and the resulting settling curves show a strong light/heavy symmetry, for the simulation ${\mathrm{C}}_{50}^{12}$ the situation is different. Due to the centrifugal effect, heavy particles are pushed away from flow vortices and light particles are pushed towards them. This results in faster settling for heavy particles and an extreme flattening of the settling curves of light particles.

Figure 15

Similar to Figs. 5 and 6, only here for the simulation set with the highest Rayleigh number: ${\mathrm{C}}_{50}^{12}$. The temperature scale is clipped, with the respective range being depicted in Fig. 16 below. The snapshot is taken at $t=0.4/\left(0.3u_{\mathrm{rms}}\right)$, i.e., at the time for which the terminal distance is 0.4 for the settling curve with $v_{\mathrm{t}}/u_{\mathrm{rms}}=0.3$ (as in Fig. 5). See Supplemental Material, Video S3 [61].

Figure 16

Demonstration of the centrifugal effect. When the convective vigor is high, particles from the transitional and bilinear range, $0.02<|v_{\mathrm{t}}|/u_{\mathrm{rms}}<2.0$, have a tendency to move towards strong flow vortices. Here we show the increased concentration of light particles with $v_{\mathrm{t}}/u_{\mathrm{rms}}=-0.3\pm 0.02$ near long-lived vortices that have developed in the simulation set ${\mathrm{C}}_{50}^{12}$. See Supplemental Material, Video S4 [61].

Figure 17

$|v_{\mathrm{t}}|t_{\mathrm{c}}/(1-{\mathcal{P}}_{\mathrm{e}})$ ratio for the simulation sets ${\mathrm{B}}_{10}^{10}, {\mathrm{C}}_{10}^{10}$, and ${\mathrm{C}}_{50}^{12}$. Light particles (circles) experience slower settling than the heavy ones, as long as the centrifugal effect comes into play (sets labeled as C). The separation between the ratios for light and heavy particles increases with the Rayleigh number (cf. the orange and green symbols).

Figure 18

Selected trajectories of particles with $v_{\mathrm{t}}/u_{\mathrm{rms}}\approx 0.4$ that have a long residence time. We show one trajectory for each of the simulation sets ${\mathrm{C}}_{50}^8, {\mathrm{C}}_{50}^{10}$, and ${\mathrm{C}}_{50}^{12}$ (blue, orange, and green lines, respectively). For clarity of the figure, the trajectories are shifted horizontally in order to avoid sidewall crossings. Sampling of the plotted trajectories is coarser than the actual numerical time step. Arrows indicate initial positions of the particles and stars mark settling events.

Figure 19

(a) Decadic logarithm of the residence time of crystals with $r_{\mathrm{c}}\in [0.1,100]$ mm and ${\rho }_{\mathrm{p}}/{\rho }_{\mathrm{f}}\in [0,2]$, suspended in a global magma ocean with parameters given in Table 1, for which we estimate $u_{\mathrm{rms}}$ to be [2.9,5.2]. The black triangles correspond to the black triangles in Fig. 10, and circles are the respective light particle types. We show ${log}_{10}$ of the value of $t_{\mathrm{t}}F$, where $t_{\mathrm{t}}$ is the terminal time $H/\left(v_{\mathrm{t}}{u}^{*}\right)$ and $F$ accounts for the characteristic settling behavior of the various particle types [see Eq. (23) and Table 2, where the parameters $\mu ,\sigma ,$ and $\lambda$ correspond to the set ${\mathrm{C}}_{50}^{12}$, i.e., to the simulation with the highest Ra]. (b) Same as above, only this time $u_{\mathrm{rms}}=0.5$ instead of 4.0 and the size of the system is $H=2.9$ km instead of $H=2890$ km. The particle set now spreads across the slow belt. We show only the heavy particles from the set (the position of particle types and the red lines are symmetrical with respect to 1 on the $x$ axis).