Linear analysis of flow mode transition triggered by finite-sized particles in Rayleigh-Bénard convection

An infinitely extended laminar system of RB convection along the $z$ direction, where the gravitational acceleration $\bm{g}$ acts in the $-\bm{e_{y}}$ direction, is considered in this study, as shown in Fig. 1(a). The cross-section of the system is a closed square cell of length $l$, with uniformly suspended spherical particles of diameter $d_{p}$ in both $x$ and $y$ direction, as shown in Fig. 1(b). The total particle number in the square cell is $N^{2}$. The lateral walls of the system are thermally insulated, and the temperature difference $\Delta T=T_{h}-T_{c}$ between the top and bottom plates is constant. Throughout this study, the following properties of the fluid and particles are regarded as constant: density $\rho$, kinematic viscosity $\nu$, thermal conductivity $\lambda$, heat diffusivity $\kappa$, volumetric thermal expansion coefficient $\beta$, and specific heat $c$. The subscripts $f$ and $p$ indicate the fluid and particle phases, respectively.

Neglecting the impact of particle dynamics on the flow motion, the nondimensional equation for the time evolution of vorticity can be given as:

$$\frac{\partial}{\partial t^{*}}\left(\nabla^{*}\times\bm{u}^{*}\right)+\nabla^{*}\times\left[(\bm{u}^{*}\cdot\nabla^{*})\bm{u}^{*}\right]=\sqrt{\frac{{\rm Pr}}{{\rm Ra}}}\nabla^{*}\times(\nabla^{*2}\bm{u}^{*})+\nabla^{*}\times(T_{f}^{*}\bm{e_{y}}),$$

(1)

where ${\rm{Pr}}=\frac{\nu_{f}}{\kappa_{f}}$ is the Prandtl number, ${\rm{Ra}}=\frac{\beta g\Delta Tl^{3}}{\nu_{f}\kappa_{f}}$ is the Rayleigh number, $\bm{u}$ is the fluid velocity, and $T_{f}$ is the fluid temperature. In the above, $\Delta T$ is the characteristic temperature, and $t_{r1}=l/U$ is set as the reference time, where $U=\sqrt{g\beta\Delta Tl}$. A coherent convective flow with negligible radial velocity is assumed to rotate around the domain center of the closed container in the cross-section and infinitely extends in the $z$ direction. Focusing on the $x\mbox{-}y$ plane, the continuity equation can be simplified to:

$$\frac{\partial u_{\theta}^{*}}{\partial\theta}=0,$$

(2)

where $\theta$ is the angular coordinate. A circular control volume of $S(r)$ bounded by $\Gamma(r)$ is supposed to cover the convective flow, as shown in Fig. 1(a), where $r$ is the radial coordinate. Based on $S$, an oriented surface with the surface element vector of $\bm{e}_{z}dS$ can be defined. This oriented surface is enclosed by an oriented curve, whose unit vector is $\bm{e}_{\theta}$, as shown in Fig. 2(b). Integrating Eq. (1) over this oriented surface gives:

$$\frac{\partial}{\partial t^{*}}\langle u_{\theta}^{*}\rangle_{\Gamma}=\sqrt{\frac{{\rm{Pr}}}{{\rm{Ra}}}}\left(\frac{\partial^{2}}{\partial r^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}+\frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\langle u_{\theta}^{*}\rangle_{\Gamma}\right)+\frac{1}{2}\left(\langle T_{f}^{*}\rangle_{\Gamma R}-\langle T_{f}^{*}\rangle_{\Gamma L}\right),$$

(3)

where $\langle\cdot\cdot\cdot\rangle_{\Gamma}$ denotes the curve average. In the above, $\langle u_{\theta}\rangle_{\Gamma}$ is the curve-averaged flow velocity, which comes from:

$$\int_{s}(\nabla\times\bm{u})\cdot\bm{e_{z}}dS=\int_{\Gamma}\bm{u}\cdot\bm{e_{\theta}}d\Gamma=\int_{\Gamma}u_{\theta}d\Gamma=\Gamma\langle u_{\theta}\rangle_{\Gamma}.$$

(4)

The convective term is simplified to zero with the continuity equation:

$$\int_{s}\nabla\times[(\bm{u}\cdot\nabla)\bm{u}]\cdot\bm{e_{z}}dS=\int_{0}^{2\pi}r\left(u_{r}\frac{\partial u_{\theta}}{\partial r}+\frac{u_{\theta}}{r}\frac{\partial u_{\theta}}{\partial\theta}\right)d\theta=0,$$

(5)

where $u_{r}$ is the radial velocity. The first term on the right-hand side of Eq. (3) comes from the viscosity term:

$$\int_{s}\nabla\times(\nabla^{2}\bm{u})\cdot\bm{e_{z}}dS=\int_{0}^{2\pi}r\left(\frac{\partial^{2}u_{\theta}}{\partial r^{2}}+\frac{1}{r}\frac{\partial u_{\theta}}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}u_{\theta}}{\partial\theta^{2}}\right)d\theta=\Gamma\left(\frac{\partial^{2}}{\partial r^{2}}\langle u_{\theta}\rangle_{\Gamma}+\frac{1}{r}\frac{\partial}{\partial r}\langle u_{\theta}\rangle_{\Gamma}\right),$$

(6)

which describes the momentum diffusion along the radial direction. The last term on the right-hand side of Eq. (3) comes from the buoyancy term:

$$\int_{s}\nabla\times(T_{f}\bm{e_{y}})\cdot\bm{e_{z}}dS=\int_{s}\frac{\partial T_{f}}{\partial x}\bm{e_{z}}\cdot\bm{e_{z}}dS=\frac{\Gamma}{2}\left(\langle T_{f}\rangle_{\Gamma R}-\langle T_{f}\rangle_{\Gamma L}\right),$$

(7)

where $S=S_{L}+S_{R}$ and $\Gamma=\Gamma_{L}+\Gamma_{R}$, as shown in Fig. 2(c). Eq. (3) shows that the imbalance of averaged fluid temperature between the left and right sides of the system is the triggering factor to invoke the bulk flow from rest to rotation.

The thermal impact of particles on flow evolution through the inter-phase heat exchange process is focused on in this study. Assuming a linear condition where the convection is suppressed and the symmetric particle array is kept without particle-particle contact during the flow motion, the evolution of fluid temperature located on the right side of the system depends mainly on the fluid-particle heat exchange:

$$\left(\rho_{f}c_{pf}\frac{\Gamma}{2}\right)\frac{d}{dt}\langle T_{f}\rangle_{\Gamma R}=-\sum_{i=1}^{N_{R}}\rho_{p}c_{pp}\frac{d}{dt}\langle T_{pi}\rangle_{vp}\Gamma_{Ri},$$

(8)

where $N_{R}$ is the total particle number distributed on $\Gamma_{R}$, $v_{p}$ is the volume of every single particle, $\langle T_{pi}\rangle_{vp}$ is the volume-averaged temperature of the ith particle on $\Gamma_{R}$, and $\Gamma_{Ri}$ is the corresponding fraction of $\Gamma_{R}$ occupied by the ith particle, as shown in Fig. 3(a). Considering the flow motion is tiny with suppressed convection, we assume the ith particle is always surrounded by the same fluid element $f_{i}$ of volume $v_{f}$ and exchanges heat only with this $f_{i}$, as shown in Fig. 3(b), where $v_{f}=\left(\frac{l}{N}\right)^{3}-v_{p}$. Therefore, the temperature evolution of the ith particle can be written as:

$$m_{p}c_{pp}\frac{d}{dt}\langle T_{pi}\rangle_{vp}=\pi d_{p}^{2}h_{p}\left[\langle T_{fi}\rangle_{vf}-\langle T_{pi}\rangle_{vp}\right],$$

(9)

where $m_{p}$ is the particle mass, $h_{p}$ is the particle heat transfer coefficient, and $\langle T_{fi}\rangle_{vf}$ is the volume-averaged temperature of the fluid around the ith particle. Combining with Eq. (8), we have:

$$\frac{\Gamma}{2}\frac{d}{dt}\langle T_{f}\rangle_{\Gamma R}=\frac{1}{H\tau_{th}}\sum_{i=1}^{N_{R}}\left[\langle T_{pi}\rangle_{vp}-\langle T_{fi}\rangle_{vf}\right]\Gamma_{Ri},$$

(10)

$$H=\frac{\rho_{f}c_{pf}}{\rho_{p}c_{pp}},$$

(11)

$$\tau_{th}=\frac{m_{p}c_{pp}}{\pi d_{p}^{2}h_{p}}=\frac{1}{6{\rm{Nu_{p}}}}\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}^{2}}{\kappa_{p}},$$

(12)

where the particle Nusselt number ${\rm{Nu_{p}}}$ comes from the dimensionless process of $h_{p}$. The empirical correlation for ${\rm Nu_{p}}$ in natural convection systems with immersed spheres can be employed [17, 18]. Similarly, the temperature evolution of $\langle T_{f}\rangle_{\Gamma L}$ can be given as:

$$\frac{\Gamma}{2}\frac{d}{dt}\langle T_{f}\rangle_{\Gamma L}=\frac{1}{H\tau_{th}}\sum_{j=1}^{N_{L}}\left[\langle T_{pj}\rangle_{vp}-\langle T_{fj}\rangle_{vf}\right]\Gamma_{Lj},$$

(13)

where $N_{L}$ is the total particle number distributed on $\Gamma_{L}$, $\langle T_{pj}\rangle_{vp}$ is the volume-averaged temperature of the jth particle on $\Gamma_{L}$, and $\Gamma_{Lj}$ is the corresponding fraction of $\Gamma_{L}$ occupied by the jth particle. Using Eq. (3), (10), and (13), we have:

	$\displaystyle\frac{\partial^{2}}{\partial t^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}$	$\displaystyle=\sqrt{\frac{{\rm{Pr}}}{{\rm{Ra}}}}\frac{\partial}{\partial t^{*}}\left(\frac{\partial^{2}}{\partial r^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}+\frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\langle u_{\theta}^{*}\rangle_{\Gamma}\right)$		(14)
		$\displaystyle+\frac{1}{\Gamma}\frac{1}{HSt_{th}}\left[\sum_{i=1}^{N_{R}}\left(\langle T_{pi}^{*}\rangle_{vp}-\langle T_{fi}^{*}\rangle_{vf}\right)\Gamma_{Ri}-\sum_{j=1}^{N_{L}}\left(\langle T_{pj}^{*}\rangle_{vp}-\langle T_{fj}^{*}\rangle_{vf}\right)\Gamma_{Lj}\right],$

where $St_{th}=\tau_{th}/t_{r1}$ is the non-dimensional time of thermal relaxation. Therefore, the flow velocity can be regarded as being updated every time interval $\Delta t^{*}=HSt_{th}$ by the inter-phase heat exchange process, which happens at the particle scale.

In this study of linear analysis, the time scale of $t_{r2}=d_{p}/U$ is defined in the particle scale. The mutually updated process of heat transfer and flow motion can be illustrated in Fig. 4, where $\Delta t=t_{r2}\Delta t^{*}$.

• •

  $t=0$: Particles and the fluid rest in thermal equilibrium.

• •

  $0<t<\Delta t$: Fluid moves with an initial perturbation velocity $\langle u_{\theta}^{*}\rangle_{\Gamma}$, resulting in a displacement of $l_{A}\equiv\langle u_{\theta}^{*}\rangle_{\Gamma}U\Delta t$ away from the thermal equilibrium position.

• •

  $t=\Delta t$ : For every particle on the right system, $\langle T_{pi}\rangle_{vp}-\langle T_{fi}\rangle_{vf}\equiv-\frac{\Delta T}{l}l_{A}=-T_{A}$. For every particle on the left system, $\langle T_{pj}\rangle_{vp}-\langle T_{fj}\rangle_{vf}\equiv\frac{\Delta T}{l}l_{A}=T_{A}$.

• •

  $\Delta t<t<2\Delta t$: Heat exchange between every particle and its surrounding fluid begins.

• •

  $t=2\Delta t$: $\langle T_{pi}\rangle_{vp}-\langle T_{fi}\rangle_{vf}$ and $\langle T_{pj}\rangle_{vp}-\langle T_{fj}\rangle_{vf}$ are updated as the value at $t=2\Delta t$.

• •

  $2\Delta t<t<3\Delta t$: Fluid moves with an updated velocity that satisfies Eq. (14).

Considering the symmetric particle arrangement on $\Gamma$ is kept, Eq. (14) can be rewritten as:

$$\frac{\partial^{2}}{\partial t^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}=\frac{\Gamma_{p}}{\Gamma}\frac{1}{HSt_{th}}\left(\langle T_{pi}^{*}\rangle_{vp}-\langle T_{fi}^{*}\rangle_{vf}\right),$$

(15)

where $\Gamma_{p}$ refers to the fraction of $\Gamma$ that occupied by the total particles distributed on $\Gamma$. Taking $t=\Delta t$ as the base state, it is suggested from the above process that the flow evolves in a way affected by the heat exchange process, where three time scales are involved:

$$\Delta t=\frac{\rho_{f}c_{pf}}{\rho_{p}c_{pp}}\frac{1}{6{\rm{Nu_{p}}}}\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}^{2}}{\kappa_{p}}\frac{d_{p}}{l},$$

(16)

$$\tau_{p}=\frac{d_{p}^{2}}{4\kappa_{p}},$$

(17)

$$\tau_{f}=\frac{d_{f}^{2}}{4\kappa_{f}},\ d_{f}=\left(\frac{6}{\pi}v_{f}\right)^{1/3}.$$

(18)

In this study, $d_{p}/l\sim 10^{-2}$ is assumed and $\tau_{f}>\Delta t$ is indicated. Therefore, the $\Delta t$ time interval is always insufficient for heat to update through the fluid region of $d_{f}$. Considering our interest in high-thermal-conductivity particles, a precondition of $\tau_{f}\gg\tau_{p}$ is further assumed, and three regimes can be defined:

• •

  I. $\tau_{p}>\Delta t$ and $\tau_{f}>\Delta t$, which means $\frac{\rho_{f}c_{pf}}{\rho_{p}c_{pp}}\frac{\lambda_{p}}{\lambda_{f}}<P$;

• •

  II. $\tau_{p}<\Delta t<\tau_{f}$, which means $\frac{\rho_{f}c_{pf}}{\rho_{p}c_{pp}}\frac{\lambda_{p}}{\lambda_{f}}>P$;

• •

  III. $\tau_{p}\sim\Delta t$ and $\tau_{f}>\Delta t$, which means $\frac{\rho_{f}c_{pf}}{\rho_{p}c_{pp}}\frac{\lambda_{p}}{\lambda_{f}}\sim P$.

In the above,

$$\frac{\tau_{p}}{\Delta t}=\frac{\rho_{p}c_{pp}}{\rho_{f}c_{pf}}\frac{\lambda_{f}}{\lambda_{p}}\frac{l}{d_{p}}\frac{3{\rm{Nu_{p}}}}{2}=\frac{\rho_{p}c_{pp}}{\rho_{f}c_{pf}}\frac{\lambda_{f}}{\lambda_{p}}P.$$

Heat exchange methods vary in different regimes, resulting in different flow modes.

In Regime I, the thermal distance traversed within the particle region during the $\Delta t$ time interval is:

$$l_{p}=\sqrt{\kappa_{p}\Delta t}<R,$$

where $R=0.5d_{p}$. Describing the thermal base state as:

$$\langle T_{pi}\rangle_{vp,t=\Delta t}-\langle T_{fi}\rangle_{vf,t=\Delta t}=T_{p0}-T_{f0}=-T_{A},$$

(19)

thermal state at the end of the first heat exchange process can be written as:

$$\langle T_{pi}\rangle_{vp,t=2\Delta t}-\langle T_{fi}\rangle_{vf,t=2\Delta t}=T_{p0}+\Delta T_{p}-T_{f0}-\Delta T_{f},$$

(20)

where $\Delta T_{p}$ refers to the volume-averaged temperature change of the single particle located on $\Gamma_{R}$ during $\Delta t<t<2\Delta t$, and $\Delta T_{f}$ refers to the volume-averaged temperature change of the corresponding fluid element. Under the initial inter-phase temperature difference, heat is transferred inside the particle region along the radial direction, which can be simplified as a one-dimensional heat conduction problem. As a result, $\Delta T_{p}$ can be written as:

$$\begin{split}\Delta T_{p}&=\frac{1}{2R}\int_{-R}^{R}d\xi\int_{-l_{p}}^{l_{p}}T_{A}\frac{{\rm{exp}}\left[-\frac{(\xi-\xi_{0})^{2}}{4\kappa_{p}\Delta t}\right]}{\sqrt{4\pi\kappa_{p}\Delta t}}d\xi_{0}\\ &=\frac{T_{A}}{2R\sqrt{4\pi\kappa_{p}\Delta t}}\frac{l_{p}}{R}\int_{-R}^{R}d\xi\int_{-R}^{R}{\rm{exp}}\left[-\frac{(\xi-\xi_{0})^{2}}{4\kappa_{p}\Delta t}\right]d\xi_{0},\end{split}$$

(21)

where $\xi$ is the integration variable referred to the existing region of every single particle on the $x\mbox{-}y$ plane, as shown in Fig. 3(b), and $\xi_{0}$ is the integration variable referred to $l_{p}$. The dual integral part in Eq. (21) can be rewritten with a transform of the integration variable:

$$\begin{split}\int_{-R}^{R}d\xi\int_{-R}^{R}{\rm{exp}}\left[-\frac{(\xi-\xi_{0})^{2}}{4\kappa_{p}\Delta t}\right]d\xi_{0}&=\int_{-R}^{R}d\xi\int_{\xi+R}^{\xi-R}{\rm{exp}}\left[-\frac{\eta^{2}}{4\kappa_{p}\Delta t}\right]d\eta\\ &=\int_{0}^{R}d\xi\left(\int_{\xi+R}^{-(\xi+R)}d\eta-\int_{-(\xi-R)}^{\xi-R}d\eta\right){\rm{exp}}\left[-\frac{\eta^{2}}{4\kappa_{p}\Delta t}\right]\\ &\approx\sqrt{4\pi\kappa_{p}\Delta t}\left(R-\int_{0}^{R}\sqrt{1-{\rm{exp}}\left[-\frac{(R-\xi)^{2}}{4\kappa_{p}\Delta t}\right]}d\xi\right),\end{split}$$

(22)

where $\eta=\xi-\xi_{0}$ . Combining with Eq. (21), we have:

$$\Delta T_{p}=\frac{T_{A}}{2R}\frac{l_{p}}{R}\left(R-\int_{0}^{R}\sqrt{1-{\rm{exp}}\left[-\frac{(R-\xi)^{2}}{4\kappa_{p}\Delta t}\right]}d\xi\right)=\frac{T_{A}}{2R}\frac{l_{p}}{R}S_{ex},$$

(23)

where $S_{ex}$ represents the integral area enclosed by $\xi=R$, $y=1$, and $y=\sqrt{1-{\rm{exp}}\left[-\frac{(R-\xi)^{2}}{4\kappa_{p}\Delta t}\right]}$, the general plot of which is shown as the red curve in Fig. 5(a). The linear part of the red curve can be similarized to fit $y_{1}=C_{1}\frac{R-\xi}{\sqrt{4\kappa_{p}\Delta t}}$, leading to an approximation of $S_{ex}$:

$$S_{ex}=\frac{1}{2}\times\frac{\sqrt{4\kappa_{p}\Delta t}}{C_{1}}\times 1=\frac{\sqrt{\kappa_{p}\Delta t}}{C_{1}},$$

(24)

where $C_{1}$ is a constant coefficient. Combining with Eq. (23), we have:

$$\Delta T_{p}=T_{A}\frac{l_{p}\sqrt{\kappa_{p}\Delta t}}{2C_{1}R^{2}}\sim T_{A}\frac{\kappa_{p}}{\kappa_{f}}\frac{d_{p}}{l},$$

(25)

$$\Delta T_{f}=-\frac{1}{HE}\Delta T_{p}\sim-T_{A}\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}}{l},$$

(26)

where $E=v_{f}/v_{p}$. Combining with Eq. (20), we have:

$$\langle T_{pi}\rangle_{vp,t=2\Delta t}-\langle T_{fi}\rangle_{vf,t=2\Delta t}\sim-T_{A}\left(1-\frac{\kappa_{p}}{\kappa_{f}}\frac{d_{p}}{l}-\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}}{l}\right).$$

(27)

Using Eq. (15), the flow motion during $2\Delta t<t<3\Delta t$ obeys the following relationship:

$$\frac{\partial^{2}}{\partial t^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}=-\frac{\Gamma_{p}}{\Gamma}\frac{d_{p}}{l}\left(1-\frac{\kappa_{p}}{\kappa_{f}}\frac{d_{p}}{l}-\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}}{l}\right)\langle u_{\theta}^{*}\rangle_{\Gamma},$$

(28)

which indicates a flow mode of simple harmonic (SH) motion, oscillating around the thermal equilibrium position.

In Regime II, the heat transfer process through the particle region should be finished for at least one round during the time interval of $\Delta t$, which means:

$$\Delta T_{p}=\frac{1}{2R}\int_{-R}^{R}d\xi\int_{-R}^{R}T_{A}\frac{{\rm{exp}}\left[-\frac{(\xi-\xi_{0})^{2}}{4\kappa_{p}\Delta t}\right]}{\sqrt{4\pi\kappa_{p}\Delta t}}d\xi_{0}=\frac{T_{A}}{2R}S_{ex},$$

(29)

$$S_{ex}=R\times 1-\frac{1}{2}\times R\times\frac{RC_{2}}{\sqrt{4\kappa_{p}\Delta t}}=R-\frac{R^{2}C_{2}}{4\sqrt{\kappa_{p}\Delta t}}.$$

(30)

In the above, $S_{ex}$ is the trapezoidal area enclosed by $\xi=R$, $\xi=0$, $y=1$, and $y=\sqrt{1-{\rm{exp}}\left[-\frac{(R-\xi)^{2}}{4\kappa_{p}\Delta t}\right]}$, the linear part of which is fitted to $y_{2}=C_{2}\frac{R-\xi}{\sqrt{4\kappa_{p}\Delta t}}$, as shown in Fig. 5(b). Therefore, the temperature change during the time interval can be written as:

$$\Delta T_{p}=T_{A}\left(\frac{1}{2}-\frac{RC_{2}}{8\sqrt{\kappa_{p}\Delta t}}\right)\sim T_{A}\left(\frac{1}{2}-\sqrt{\frac{\kappa_{f}}{\kappa_{p}}\frac{l}{d_{p}}}\right),$$

(31)

$$\Delta T_{f}=-\frac{1}{HE}\Delta T_{p}\sim-T_{A}\left(\frac{1}{2HE}-\frac{\lambda_{p}}{\lambda_{f}}\frac{\kappa_{f}}{\kappa_{p}}\sqrt{\frac{\kappa_{f}}{\kappa_{p}}\frac{l}{d_{p}}}\right).$$

(32)

Combining with Eq. (15) and Eq. (20), the flow motion during $2\Delta t<t<3\Delta t$ obeys the following relationship:

$$\frac{\partial^{2}}{\partial t^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}=-\frac{\Gamma_{p}}{\Gamma}\frac{d_{p}}{l}\left(\frac{1}{2}-\frac{1}{2HE}+\sqrt{\frac{\kappa_{f}}{\kappa_{p}}\frac{l}{d_{p}}}+\frac{\lambda_{p}}{\lambda_{f}}\frac{\kappa_{f}}{\kappa_{p}}\sqrt{\frac{\kappa_{f}}{\kappa_{p}}\frac{l}{d_{p}}}\right)\langle u_{\theta}^{*}\rangle_{\Gamma},$$

(33)

which also indicates a flow mode of SH motion around the thermal equilibrium position.

After the above discussion, the flow pattern in Regime III, where $\Delta t$ is about the time to complete one round of heat transfer through the particle region, becomes easier to understand. As shown in Fig. 5(c), the averaged temperature change can be given as:

$$\Delta T_{p}=\frac{T_{A}}{2R}S_{ex}=\frac{T_{A}}{4},$$

(34)

$$\Delta T_{f}=-\frac{1}{HE}\Delta T_{p}=-\frac{T_{A}}{4HE}.$$

(35)

As a result, flow motion during $2\Delta t<t<3\Delta t$ obeys the following relationship:

$$\frac{\partial^{2}}{\partial t^{*2}}\langle u_{\theta}^{*}\rangle_{\Gamma}=-\frac{\Gamma_{p}}{\Gamma}\frac{d_{p}}{l}\left(\frac{3}{4}-\frac{1}{4HE}\right)\langle u_{\theta}^{*}\rangle_{\Gamma},$$

(36)

which indicates the flow period has no obvious connection with the thermal conductivity ratio.

Regime III occurs only when $\sqrt{\kappa_{p}\Delta t}\approx R$, which means it is easy to shift to Regime I or Regime II. The previous discussion shows that with the increase of $\lambda_{p}/\lambda_{f}$, the oscillation period in Regime I becomes shorter while the oscillation period in Regime II becomes longer. Therefore, there should exist a SH flow mode in Regime III, which yields the minimum oscillation period.

Considering the symmetric particle array is always kept in this study, it is reasonable to assume that particles distributed on $\Gamma$ have the same velocity:

$$u_{pi,\Gamma}=u_{p,\Gamma}=u_{p\theta,\Gamma},$$

where $u_{p\theta,\Gamma}$ is the particle velocity along the angular direction. Therefore, the drag force experienced by every single particle on $\Gamma$ can be written as:

$$f_{d}=m_{p}\frac{\langle u_{\theta}\rangle_{\Gamma}-u_{p\theta,\Gamma}}{\tau_{m}},$$

(37)

where a positive $f_{d}$ indicates the drag force acts along the counterclockwise direction, and a negative $f_{d}$ indicates the drag force acts along the clockwise direction. In the above, $\tau_{m}$ is the particle relaxation time, which can be written in the general form[19]:

$$\frac{1}{\tau_{m}}=\frac{3}{4}C_{D}f(\epsilon)^{m}\frac{\rho_{f}}{\rho_{p}}\epsilon^{2}\frac{\left|\langle u_{\theta}\rangle_{\Gamma}-u_{p\theta,\Gamma}\right|}{d_{p}},$$

(38)

where $C_{D}$ is the drag coefficient for single sphere suspension proposed by Schiller and Naumann [20], $\epsilon$ is the local void fraction or porosity, and $f(\epsilon)^{m}$ is the porosity function to correct the relaxation time when the particle suspension is not dilute enough.

Treating particles distributed on $\Gamma$ as symmetric particle couples, for example, particle $A$ and $B$ in Fig. 6, the moment of force experienced by every particle couple can be written as:

$$\sum M_{AB}=(\rho_{f}v_{p}ga-\rho_{f}v_{p}ga)+(-m_{p}ga+m_{p}ga)+2f_{d}r=2f_{d}r,$$

(39)

suggesting the particle motion is mainly driven by the drag force. According to the conservation of angular momentum, we have:

$$\frac{d}{dt}L_{p,\Gamma}=\frac{N_{\Gamma}}{2}I_{p}\frac{dw_{\Gamma}}{dt}=\frac{N_{\Gamma}}{2}\sum M_{AB},$$

(40)

where $L_{p,\Gamma}$ represents the angular momentum of particles located on $\Gamma$, $N_{\Gamma}$ is the total particle number on $\Gamma$, $I_{2}$ is the moment of inertia of every particle couple, and $w_{\Gamma}$ is the angular velocity of every particle. Assuming the particle movement is tiny enough that $w_{\Gamma}\approx u_{p\theta,\Gamma}$, Eq. (43) can be rewritten as:

$$\frac{d}{dt}u_{p\theta,\Gamma}=\frac{2rm_{p}}{I_{p}\tau_{m}}\left(\langle u_{\theta}\rangle_{\Gamma}-u_{p\theta,\Gamma}\right),$$

(41)

which suggests that particles also move with a kind of regular oscillation mode driven by the flow, while a phase difference exists between the particle SH and fluid SH motions.

An oscillatory granular flow, which is driven by the simple harmonic fluid motion, is predicted from the above analysis. In this section, validation is carried out by comparing the oscillation period suggested by the current study with that from related simulation research.

As shown in Fig. 7, red dots refer to the simulation results of particle motion in the solid-dispersed Rayleigh-Bénard convection[16], where $d_{p}/l$ is set as 0.05 and $H$ is set as 1. Fitting these conditions, Regime I of this study falls into the range of $\frac{\lambda_{p}}{\lambda_{f}}<60$ , and the granular flow yields a SH motion, whose period is:

$$\tau_{I}=\left[\frac{\Gamma_{p}}{\Gamma}\frac{d_{p}}{l}\left(1-H\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}}{l}-\frac{\lambda_{p}}{\lambda_{f}}\frac{d_{p}}{l}\right)\right]^{-\frac{1}{2}}\sim(\lambda_{p}/\lambda_{f})^{-\frac{1}{2}},$$

(42)

showing agreement with the simulation study. Regime II falls in the range of $\frac{\lambda_{p}}{\lambda_{f}}>60$ , where the granular flow yields a regular oscillation period:

$$\tau_{II}=\left[\frac{\Gamma_{p}}{\Gamma}\frac{d_{p}}{l}\left(\frac{1}{2}-\frac{1}{2HE}+\sqrt{\frac{\kappa_{f}}{\kappa_{p}}\frac{l}{d_{p}}}+\frac{\lambda_{p}}{\lambda_{f}}\frac{\kappa_{f}}{\kappa_{p}}\sqrt{\frac{\kappa_{f}}{\kappa_{p}}\frac{l}{d_{p}}}\right)\right]^{-\frac{1}{2}}\sim\left(\lambda_{p}/\lambda_{f}\right)^{\frac{1}{4}}.$$

(43)

This correlation between the oscillation period and thermal conductivity ratio shows especially good agreement with the simulation study. Regime III, which is between Regime I and II, falls in the range around $\frac{\lambda_{p}}{\lambda_{f}}=60$, where a minimum oscillation period can be observed. In the simulation results, the shortest oscillation period of all does occur around $\lambda_{p}/\lambda_{f}=60$ while the oscillation period still shows a relationship with the thermal conductivity ratio.

It is also noted that the correlation in Regime I fails to cover the flow motion around the range of $1<\frac{\lambda_{p}}{\lambda_{f}}<10$, as indicated from the simulation results. The reason is supposed to be that the assumption of $\tau_{f}\gg\tau_{p}$ becomes invalid under the problem setting of $\frac{\lambda_{p}}{\lambda_{f}}<10$, and it is not scientific to simplify the particle temperature change into a Cauchy problem anymore. More heat transferred through the particle surface to the fluid element should be modeled.