2.1. Gradual density gradient

Here, we consider a two-dimensional (2-D) ($x$–$y$ plane) incompressible system where density and potential gradients are gradually increasing along the y-direction while a constant gravitational force acts in the opposite direction. We have assumed that the charge density of dusty plasma remains constant, i.e. ${\rho _c}={\rho _{c0}}$. In the laboratory dusty plasma experiments, the external electric field is used to levitate/equilibrate dust grains at a particular height against the gravitational pull of the earth. Mathematically, we obtain such a static equilibrium situation, with some assumptions such as $\boldsymbol v_{d0} =0$, i.e. no initial flow and ${\partial }/{\partial {t}} =0$, from the above (2.1), (2.2) and (2.4),

(2.8)\begin{equation} {\rho_{d0}}{g} ={-}{\rho_{c0}}{\frac{\partial\phi_{d0}} {\partial y}}. \end{equation}

A small perturbation in the various fields, e.g. density, scalar potential and dust velocity, can be written as

(2.9)\begin{gather} {\rho_d}(x,y,t)={\rho_{d0}}(y,t=0)+{\rho_{d1}}(x,y,t), \end{gather}

(2.10)\begin{gather} {\phi_d}(x,y,t)={\phi_{d0}}(y,t=0)+{{\phi_{d1}}}(x,y,t), \end{gather}

(2.11)\begin{gather} {\boldsymbol{v}_d}(x,y,t)=0+{\boldsymbol{v}_{d1}}(x,y,t), \end{gather}

respectively. Retaining only linear terms in the perturbed fields in (2.1), (2.2) and (2.4) we obtain the following equations for the linearized instability analysis:

(2.12)\begin{equation} \frac{\partial \rho_{d1} }{\partial t} + \left(\boldsymbol{v}_{d1}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\rho_{d0}=0, \end{equation}

${y}$ component

(2.13)\begin{equation} \left[1 + \tau_m \frac{\partial} {\partial t}\right] \left[{{\rho_{d0}}{\frac{\partial {v}_{d1y}} {\partial t}}}+{\rho_{d1}}g +{\rho_{c0}}{\frac{\partial \phi_{d1y}}{\partial y}}\right]=\eta {\nabla}^2 {v}_{d1y}, \end{equation}

${x}$ component

(2.14)\begin{gather} \left[1 + \tau_m \frac{\partial} {\partial t}\right] \left[{{\rho_{d0}}{\frac{\partial {v}_{d1x}} {\partial t}}} +{\rho_{c0}}{\frac{\partial \phi_{d1x}}{\partial x}}\right]=\eta {\nabla}^2 {v}_{d1x}, \end{gather}

(2.15)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol v_{d1}=0. \end{gather}

Since equilibrium fields vary along ${y}$, the above set of equations can only be Fourier analysed in time and spatial coordinate $x$. However, we assume the perturbations to be of much smaller scale compared with the equilibrium variation along $y$. So, here, we invoke the local approximation and proceed with the Fourier decomposition along y also along with $x$. For this decomposition, the perturbed quantities are considered proportional to $\exp (\textrm {i}{k_x}x+\textrm {i}{k_y}y-\textrm {i}{\omega }{t})$. This gives relations ${\partial }/{\partial {t}}={-\textrm {i}{\omega }}$, ${\partial }/{\partial {x}}= {\textrm {i}{k}_x}$ and ${\partial }/{\partial {y}}={\textrm {i}{k}_y}$. Now, replacing the time and spatial derivatives of the perturbed quantities in (2.12), (2.13), (2.14) and (2.15) with the these relations, we get

(2.16)\begin{gather} {\rho_{d1} }={-}\frac{\textrm{i}}{{\omega}}\frac{\partial \rho_{d0}}{\partial y}{v_{d1y}}, \end{gather}

(2.17)\begin{gather} {\left[1-\textrm{i}{\omega}\tau_{m}\right]} {\left[{-\textrm{i}{\omega}{\rho_{d0}}{{v}_{d1y}}} +{\rho_{d1}}g+\textrm{i}{{k}_y}{\rho_{c0}}{\phi_{d1}}\right]}={-}{\eta {{k}}^2{v}_{d1y}}, \end{gather}

(2.18)\begin{gather} {\left[1-\textrm{i}{\omega}\tau_{m}\right]}{\left[{-\textrm{i}{\omega}{\rho_{d0}}{{v}_{d1x}}}+\textrm{i}{ {k}_x}{\rho_{c0}}{\phi_{d1}}\right]} ={-}{\eta {{k}}^2 {v}_{d1x}}, \end{gather}

(2.19)\begin{gather} \textrm{i}{{k}_x}{v_{d1x}}={-}\textrm{i}{{k}_y}{v_{d1y}} \Rightarrow {v_{d1y}}={-\frac{{k}_x}{{k}_y}}{v_{d1x}}. \end{gather}

Here, ${k}=({k_x}, {k_y})$ is the wavenumber associated with the perturbed quantities. Obtaining $\phi _{d1}$, using above relation and (2.16) (2.17) (2.18) and (2.19), we get

(2.20)\begin{equation} {{\phi_{d1}}} ={-}{\frac{{g}{{k}_x}{v_{d1x}}}{{\rho_{c0}}{\omega}{{k}^2}}} \frac{\partial \rho_{d0} }{\partial y}, \end{equation}

and using the above relation (2.20) in (2.18), we get the dispersion relation as

(2.21)\begin{equation} {\left(1-\textrm{i}{\omega}\tau_{m}\right)}{\left[{\omega^2}{{k}^2}+ {\frac{g{k}^2_x}{{{\rho_{d0}}}}} \frac{\partial \rho_{d0} }{\partial y}\right]}={-\textrm{i}{\omega}{\eta^*}{{k}}^4}, \end{equation}

where kinematic viscosity $\eta ^*={\eta }/{\rho _{d0}}$, the ratio of absolute viscosity ${\eta }$ to density ${\rho _{d0}}$. A similar equation for the electron ion plasma in the context of inertial fusion has been predicted by Das & Kaw (Reference Das and Kaw2014). In the weak coupling limit, ${\omega }{\tau _m}\ll 1$, the instability growth rate shows a viscous damping due to the dust–dust collisions. In the strongly coupled limit, i.e. ${\omega }{\tau _m}\gg 1$, the above dispersion (2.21) becomes

(2.22)\begin{equation} {\omega^2}={\frac{\eta^*}{{\tau_{m}}}{{k}}^2}-{\frac{g}{{{\rho_{d0}}}}} \frac{\partial \rho_{d0} }{\partial y}{\frac{{{k^2_x}}}{{{k}}^2}}. \end{equation}

Here, it should be noted that the analytical dispersion relation for dusty plasma medium has also been obtained by Avinash & Sen (Reference Avinash and Sen2015). However, we had presented this result earlier in 2014 in a conference (Dharodi & Das Reference Dharodi and Das2014).

2.2. Sharp interface

For a sharp density interface, the dispersion relation is obtained by replacing the gravity driven term $-({g/{\rho _{d0}}})({\partial {\rho _{d0}}}/{\partial {y}})({{{k^2_x}}/{{{k}}^2}})$ in (2.22) by $-g{{k}}{{A_T}}$ (Avinash & Sen Reference Avinash and Sen2015),

(2.23)\begin{equation} {\omega^2}={\frac{\eta^*}{{\tau_{m}}}{{{k}}^2}}-g{{k}}{{A_T}}, \end{equation}

where ${{A_T}}={(\rho _{d02}-\rho _{d01})}/{(\rho _{d01}+\rho _{d02})}$ is the Atwood number that represents the non-dimensional density difference between two adjacent fluids with a common interface. Here $\rho _{d02}$ ($\rho _{d01}$) is the density of the heavier (lighter) fluid.

In the right-hand side of (2.22) and (2.23), the second term assists the growth of GD instabilities which grow at the maximum rate in the inviscid limit ($\eta =0$) of the medium while the first shear wave term suppresses it. As our interest is in a condition for which a medium favours the GD instabilities, then for all the numerical simulations we have chosen the parametric values of all variables in such a way that the second term dominates over the first term, with a constant value of gravitational acceleration $g(= 10)$. While this condition is ensured, one has the freedom to choose different values of parameters, including the value of $g$, since the phenomenon of GD instability will take place but may be at different rates. Equations (2.22) and (2.23) with the assumption $g=0$ become (2.7), which tell us that for the fixed ratio of ${\eta }/{\tau _m}$, the phase velocity of TS waves is inversely proportional to the square root of the density of the medium. Recently, Dharodi (Reference Dharodi2020) has shown numerically that in a two-fluid medium, TS waves travel slower in the denser side than in the lighter side. This difference in speeds induced a net flow of the medium towards the lighter side. In § 3.2, this numerical result will help us to understand the role of TS waves in initiating the RT instability for VE fluids.

3.1. Simulation methodology

For the purpose of numerical simulation the generalized momentum (2.2) has been expressed as the following set of two coupled convective equations:

(3.1)\begin{gather} {{\rho_d}\left(\frac{\partial{\boldsymbol{v}_d}}{\partial {t}}+{\boldsymbol{v}_d}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol{v}_d}\right)}+{\rho_d}\boldsymbol{g} +{\rho_c}\boldsymbol{\nabla} \phi_{d}={\boldsymbol{\psi}}, \end{gather}

(3.2)\begin{gather} \frac{\partial {\boldsymbol{\psi}}} {\partial t}+\boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol{\psi}}={\frac{\eta}{\tau_m}}{{\nabla}^2}{\boldsymbol{v}_d }-{\frac{\boldsymbol{\psi}}{\tau_m}}. \end{gather}

For 2-D studies all the variables are functions of $x$ and $y$ only, i.e. ${\boldsymbol {\psi }}(x,y)$, ${\boldsymbol v_d}(x,y)$, ${\rho _d}(x,y)$ and $\boldsymbol {g}(y)$. From (3.1) it is clear that the quantity ${\boldsymbol {\psi }}(x,y)$ is the strain created in the elastic medium by the time-varying velocity fields. We rewrite (3.1) using the equilibrium condition (2.8) in terms of the perturbed density relation (2.9):

(3.3)\begin{equation} \frac{\partial{\boldsymbol{v}_d}}{\partial {t}}+{\boldsymbol{v}_d}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol{v}_d}+{\frac{\rho_{d1}}{\rho_d}}\boldsymbol{g} +{\frac{\rho_{c0}}{\rho_d}}{\boldsymbol{\nabla}{\phi_{d1}}}=\frac{\boldsymbol{\psi}}{\rho_d}. \end{equation}

Taking the curl of (3.3) and using the Boussinesq approximation, we get

(3.4)\begin{equation} \frac{\partial{\xi_{z}}} {\partial t}+\left(\boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla}\right) {\xi_{z}}={\frac{1}{\rho_{d0}}}{{\boldsymbol{\nabla}}\times{\rho_{d1}}{\boldsymbol{g}}}+{\boldsymbol{\nabla}}\times{\frac{\boldsymbol{\psi}}{\rho_d}}. \end{equation}

Here, ${\xi _{z}}(x,y)=[{\boldsymbol {\nabla }} \times {\boldsymbol v_d}(x,y) ]_z$ is the $z$ (and the only finite) component of vorticity for the 2-D case considered here. It is normalized by the dust plasma frequency. The acceleration $\boldsymbol g$ is applied opposite to the fluid density gradient $-g\hat {y}$. The coupled set of model equations for numerical simulations are the following:

(3.5)\begin{gather} \frac{\partial \rho_d }{\partial t} + \left(\boldsymbol{v}_d\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\rho_d=0, \end{gather}

(3.6)\begin{gather} \frac{\partial {\boldsymbol{\psi}}} {\partial t}+\left(\boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla}\right) {\boldsymbol{\psi}}={\frac{\eta}{\tau_m}}{{\nabla}^2}{\boldsymbol{v}_d }-{\frac{\boldsymbol{\psi}}{\tau_m}}, \end{gather}

(3.7)\begin{gather} \frac{\partial{\xi}_z} {\partial t}+\left(\boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla}\right) {{\xi}_z}={-}{\frac{g}{\rho_{d0}}}{\frac{\partial{\rho_{d1}}} {\partial x}} +{\frac{\partial}{\partial x}}\left({\frac{\psi_{y}}{\rho_d}}\right) -{\frac{\partial}{\partial y}}\left({\frac{\psi_{x}}{\rho_d}}\right). \end{gather}

For the numerical simulations, we have used the LCPFCT method (Boris et al. Reference Boris, Landsberg, Oran and Gardner1993) to evolve the coupled set of (3.5), (3.6) and (3.7). This method is based on a finite difference scheme associated with the flux-corrected algorithm. The velocity at each time step is updated by using Poisson's equation ${{\nabla }^2}{\boldsymbol {v}_d}=-{{\boldsymbol {\nabla }}}\times {\boldsymbol {\xi }}$, which has been solved using the FISPACK package (Swarztrauber, Sweet & Adams Reference Swarztrauber, Sweet and Adams1999). Here, it should be noted that the equations for simulations have been derived under the Boussinesq approximation, which was not assumed in the above analytical derived dispersion relations ((2.22) and (2.23)). So, our numerical simulations are used to show the same trend of GD instabilities growth as expected analytically instead of quantitative matching. Throughout simulation studies, boundary conditions are periodic in the horizontal direction ($x$-axis) while non-periodic along the vertical ($y$-axis) direction where the effects of perturbed quantities die out before hitting the boundary of the simulation box. Moreover, we have stopped all simulations for both directions before hitting the boundaries. In each case, the grid convergence study has been carried out first to make sure of grid independence of the numerical results and to avoid the simulation resolution issue.

In the interest of observing the pure viscous damping effect on the growth of instabilities, we have also simulated the cases of pure viscous HD fluids with the value of $\tau _m=0$. For such cases (3.6) has a singularity. To avoid such a singularity, we put ${\tau _m}=0$ into the generalized momentum (2.2) and we obtain the standard HD equation (2.5). Take the curl of (2.5), under the same above considered conditions and assumptions, we get the HD vorticity equation,

(3.8)\begin{equation} \frac{\partial{\xi}_z} {\partial t}+\left(\boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla}\right) {{\xi}_z}={-}{\frac{g}{\rho_{d0}}}{\frac{\partial{\rho_{d1}}} {\partial x}}+{\eta}{{\nabla}^2}{{\xi}_z}. \end{equation}

Thus, numerically, for pure HD cases ($\tau _m=0$) we solve the set of (3.5) and (3.8) where dust fluid velocity at each time step is updated by using ${{\nabla }^2}{\boldsymbol {v}_d}=-{{\boldsymbol {\nabla }}}\times {\boldsymbol {\xi }}$. Since we have simulated the normalized equations, so in all figures the numbers in colour bars correspond to the normalized value of that plotted physical quantity.

3.2. RT instability

We consider two types of inhomogeneous density profiles, namely case A with a sharp interface between two different densities where heavier fluid $\rho _{d02}$ (upper) is supported by lighter one $\rho _{d01}$ (lower) as shown in figure 1(a), and case B with gradually increasing density gradient along the $y$-axis is given by (3.10) that is shown in figure 1(b). For both profiles, numerical simulations are performed with $512\times 512$ grid points in both the $x$ and $y$ directions on a simulation region in the $x$–$y$ plane defined by $-{{\rm \pi} } < x < {{\rm \pi} }$, $-{{\rm \pi} } < y < {{\rm \pi} }$. Thus, we have a system of the length of $lx=ly=2{{\rm \pi} }$ units. The mathematical formalism of the initial equilibrium condition can be understood through the relation (2.8), a sufficient electric field has been applied against the gravity to levitate the dust grains.

Figure 1. Initial density profiles at $t=0$ to study the RT instability. The patches along the $x$ direction at the centre of vertical axis represent the imposed sinusoidal perturbation which is given by (3.9) to hasten the instability. (a) Sharp density interface profile and (b) gradually increasing density gradient along the $y$ direction, or say against the gravity; where $\boldsymbol {H}$ stands for heavy density and $\boldsymbol {L}$ for the lighter one.

3.3. Sharp interface

For case A (figure 1a), we consider a system which consists of two incompressible fluids of constant densities $\rho _{d01}=1$ for $-{{\rm \pi} }\leq y \leq 0$ (lower half) and $\rho _{d02}=2$ for $0=y \leq {{\rm \pi} }$ (upper half), i.e. the denser fluid of density $\rho _{d02}$ is positioned above the lighter fluid $\rho _{d01}$. This is a very commonly used density profile to study RT instability; a few references are Waddell et al. (Reference Waddell, Niederhaus and Jacobs2001), Bell et al. (Reference Bell, Day, Rendleman, Woosley and Zingale2004), Avinash & Sen (Reference Avinash and Sen2015), Hosseini et al. (Reference Hosseini, Turek, Möller and Palmes2017) and Liang et al. (Reference Liang, Hu, Huang and Xu2019). In the absence of gravity, this system will remain stable. To hasten the evolution of our stable system under gravity, we impose a small amplitude ${\bar {\rho }_{d1}}$ sinusoidal perturbation (Tiwari et al. Reference Tiwari, Das, Angom, Patel and Kaw2012) on the density profile at the interface ($y=0$), i.e.

(3.9)\begin{equation} \rho_{d1}={\bar{\rho}_{d1}}{\cos({{k}}_xx)}\exp({-}y^2/{\sigma^2}). \end{equation}

Here, the Gaussian function $\exp (-y^2/{\sigma ^2})$ controls the width of perturbation at the interface $y=0$ along the $y$-axis through the parameter $\sigma$, the larger the value of $\sigma$ the greater the width of the perturbation; $\sigma$ is supposed to be ${\ll }l_y/2$, half of the box size. We have chosen $\sigma =0.1$, the maximum amplitude of perturbation ${\bar {\rho }_{d1}}=0.01$, and the longest possible single mode ${{k}}_x=1$ along the $x$ direction (${{k}}_x=N({2{\rm \pi} }/{lx})$, the positive integer $N$ represents the number of the modes, here, $N=1$). The total density is $\rho _d=\rho _{d0s}+\rho _{d1}$, where the suffix $s$ corresponds to $1$ and $2$, which are labelled for the low and high density dust fluid, respectively.

The evolution of the RT instability is reproduced for the inviscid HD fluid from our code (see figure 2). The initial perturbation was given at ${{k}}_x=1$, which is observed. The heavy fluid is observed to penetrate into the lighter fluid and the lighter fluid rises above into the domain of the heavy fluid. Subsequently, the falling heavy fluid develops rolls at the edges acquiring a mushroom-like profile. These rolls grow up with time as the heavy fluid keeps falling with time. Such an evolution of RT instability is a well known result of HD fluid (Hosseini et al. Reference Hosseini, Turek, Möller and Palmes2017; Talat et al. Reference Talat, Mavrič, Hatić, Bajt and Šarler2018), thereby substantiating our numerical simulation.

Figure 2. Time evolution of RT instability for inviscid HD fluid at the sharp interface of two different densities. The heavy fluid (upper half) penetrates into the lighter fluid (lower half) and the lighter fluid rises above into the domain of heavy fluid. Subsequently, the relative flow of two adjacent fluids starts rolling which causes the formation of vortex-like structures at the interface. The final structure acquires a mushroom-like profile.

In order to notice the pure viscous damping effect of viscosity on the growth of the RT instability, we have also simulated a pure viscous HD fluid (see the density evolution in figure 3 for $\eta = 0.1$, $\tau _m=0$) by using (3.5) and (3.8). In the viscous fluid (figure 3) the penetration of heavy fluid into the lighter fluid is much slower than with inviscid fluid (figure 2) that results in no rolls at the interface even at almost double the time step of the inviscid case. Next, this viscous fluid (figure 3) has been compared with a VE fluid, shown in figure 4(a–d), which has the same viscosity $\eta =0.1$ but a finite value of $\tau _m =20$. We have discussed that in VE fluids $\eta$ in the presence of $\tau _m$ assists the TS waves which travel slower in the denser fluid than in the lighter. In our two fluid case, the lighter fluid is on the lower side while the denser fluid is on the upper, which means the difference in speeds of these waves should favour the net flow of fluid in the downward direction. This observation can be clearly seen from the comparison of two mediums where the piercing of heavy fluid into the lighter fluid is faster in a VE fluid (figure 4a–d) than a pure viscous fluid (figure 3). Hence, in the present scenario, the RT instability grows faster in a VE fluid in comparison with a viscous fluid. Moreover, this comparison shows that the instability grows phenomenologically similar in both fluids as observed in the dusty plasma experiment by Pacha et al. (Reference Pacha, Heinrich, Kim and Merlino2012).

Figure 3. Time evolution of RT instability for pure viscous HD fluid ($\eta = 0.1$; $\tau _m=0$) at the sharp interface of two different densities. Here, due to the viscous damping force, the RT instability grows at a considerably slower pace compared with the inviscid HD fluid (see figure 2). Numerically, the set of (3.5) and (3.8) has been solved to avoid the singularity in generalized momentum (2.2) due to ${\tau _m}=0$.

Figure 4. The growth of RT instability at the sharp interface of two VE fluids of different densities: (a–d) $\eta =0.1$, $\tau _m=20$; (e–h) $\eta =0.1$, $\tau _m=5$. The RT instability gets weaker with increasing ratio ${\eta }/{\tau _m}$ as expected analytically from (2.23). This can be noticed by the comparison of panels (a–d) and (e–h), enlargement in width of the neck region (as indicated by double sided arrow) and reduction in fluid flow in the downward direction at each time step.

Contour plots in figure 4 show that the growth rate of RT instability get suppressed with increasing coupling strength (ratio $\eta$/$\tau _m$) of the medium. The evolution of RT instability or the downward fluid flow under an applied gravitational force causes the appearance of free energy in the medium. This free energy release in the form of TS waves which travel with velocity $\sqrt {\eta /\tau _m}$ in the perpendicular direction to the evolution of the instability that, in turn, slow down the fluid flow. These transverse effects get stronger with an increase of the ratio $\eta /\tau _m$ that can be noticed by the comparison of figures 4(a–d) and 4(e–h), enlargement in width of the neck region (as indicated by double sided arrow) and reduction in fluid flow in the downward direction at each time step. Furthermore, the comparative evolution of density in figures 3, 4(a–d) and 4(e–h) shows that the appearance of $\tau _m$ speed up the growth of RT instability ($\eta =0.1$ is same for all the three cases). Both the above observations favour the analytical expected result from (2.23).

3.4. Gradual variation of density profile

Next, the density of the dust fluid is stratified against gravity and has a gradually changing profile as shown in figure 1(b). We assume that the density has an exponential shape in $y$ (such a profile has been discussed in chapter 19 in Goldston (Reference Goldston2020)), i.e.

(3.10)\begin{equation} \rho_{d}={\rho_{d0}}+{\bar{\rho}_{d0}}\exp({\epsilon}y)+{\rho_{d1}}. \end{equation}

Here, $\rho _{d0}$ is constant background density and ${\bar {\rho }_{d0}}$ is the amplitude of inhomogeneity. The parameter $\epsilon$ decides the ramp of inhomogeneity in background density. The values of parameters taken for the present case are $\rho _{d0}=5$, ${\bar {\rho }_{d0}}=0.1$ and $\epsilon =0.025$. The sinusoidal perturbation $\rho _{d1}$ (see (3.9) with ${{\bar {\rho }_{d1}}}=0.01$, $\sigma =0.1$ and ${{k}}_x=1$) is imposed at $y = 0$, which can be clearly seen in figure 1(b).

In figure 5, we have compared two different VE fluids with the same $\eta$, (${\eta }=0.1, {\tau _m}=20, {\eta }/{\tau _m}=0.005$ for figure 5a–d; ${\eta }=0.1, {\tau _m}=5, {\eta }/{\tau _m}=0.02$ for figure 5e–h). From the comparison of figure 5(a–d) and figure 5(e–h), it is clear that at each time step, the growth of RT instability gets weaker with increasing coupling strength, while ${\tau _m}$ speeds up it as suggested by (2.22); similar observations have been made for the sharp density gradient profiles.

Figure 5. The time evolution of RT instability for two VE fluids having gradually increasing densities against the gravity. (a–d) The VE fluid with $\eta =0.1$, $\tau _m=20$ and (e–h) VE fluid with $\eta =0.1$, $\tau _m=5$. The RT instability gets weaker with increasing coupling strength, i.e. ${\eta }/{\tau _m}$ as analytically suggested by (2.22).

3.5. BD evolution

We consider the two types of density inhomogeneity profiles, namely case A an initially static and circular bubble whose density is less than that of the surrounding quiescent HD/VE fluid as shown in figure 6(a), and case B initially static and a circular droplet whose density is higher than that of the surrounding medium as shown in figure 6(b). In order to avoid the confusion between the evolution of density and evolution of vorticity, we have discussed density evolution in terms of blob/blobs while vorticity evolution in terms of  lobe/lobes. For both cases, we have considered a system of length $lx=ly=12{{\rm \pi} }$ units with $512\times 512$ grid points in both the $x$ and $y$ directions. The total density is $\rho _d=\rho _{d0}+\rho _{d1}$, where ${\rho _{d0}}$ is the background density and the inhomogeneous Gaussian density $\rho _{d1}$ of rising bubble or falling droplet centred at $({x_{c}},{y_{c}})$, with radius ${a_{c}}$ is

(3.11)\begin{equation} \rho_{d1}={\bar{\rho}_{d1}}\exp\left(\frac{-\left({\left(x-x_{c} \right)^2+(y-y_{c})^2}\right)}{a^2_{c}}\right). \end{equation}

Here $a_c$ is equal for both blobs. The parameter ${\bar {\rho }_{d1}}$ has a negative value for the bubble while a positive value of equal magnitude for droplet. Thus, in both cases the variation in density is symmetric about an axis that is perpendicular to the simulation plane.

Figure 6. Initial densities’ profiles at time $t=0$. (a) Bubble of less density inside the heavier fluid and (b) droplet of high density inside the lighter fluid. In the colour bar, letter $\boldsymbol {{H}}$ is the heavy density region and $\boldsymbol {L}$ stands for the lighter density region.

Before running codes and reporting simulation results merely based on the interpretations of contour plots, let us discuss the equations of the model to develop some simple qualitative understanding of simulations results. For an inviscid flow, the vorticity (3.4), (3.7), (3.8) in HD limit, becomes

(3.12)\begin{equation} \frac{\partial{\xi_{z}}} {\partial t}+\left(\boldsymbol{v}_d \boldsymbol{\cdot} \boldsymbol{\nabla}\right) {\xi_{z}}={\frac{1}{\rho_{d0}}}{{\boldsymbol{\nabla}}\times{\rho_{d1}}{\boldsymbol{g}}}={-}{\frac{g}{\rho_{d0}}}{\frac{\partial{\rho_{d1}}} {\partial x}}. \end{equation}

This equation tells that as the simulation begins, the influence of gravity induces a dipolar vorticity. The dipolar vorticity, $-{({g}/{\rho _{d0}})}{({\partial {\rho _{d1}}}/{\partial x})}={({2g}/{\rho _{d0}})}{(x-x_{c})}{\rho _{d1}}$, acts as a buoyant force on droplet, while dipolar vorticity, $-{({g}/{\rho _{d0}})}{({\partial {\rho _{d1}}}/{\partial x})}=-{({2g}/{\rho _{d0}})}{(x-x_{c})}{\rho _{d1}}$ acts as a buoyant force on bubble. Both these are oppositely propagating dipolar vorticities and each have a pair of counter-rotating lobes (or a pair of unlike-sign lobes). The resultant of a counter-rotating lobe pair related to the bubble causes vertical upward motion, and vertical downward motion for the droplet (Horstmann et al. Reference Horstmann, Henningsson, Thomas and Bomphrey2014). It can be visualized through the schematic diagram of vorticity in figure 7, corresponding to the bubble/droplet density in figure 6. The curved solid arrows over the lobes represent the direction of rotation of lobes, while the net propagation is indicated by vertical dotted arrows.

Figure 7. Schematic vorticity diagram (a pair of counter-rotating lobes, the rotational direction of each lobe is shown by a curved arrow) corresponding to figure 6. (a) Dipolar vorticity due to the buoyant force on a bubble (figure 6a). The net rotation causes the motion in the vertical upward direction, indicated by the vertical upward dotted arrow. Similarly, vorticity in panel (b) with respect to figure 6(b) causes the net downward flow for the droplet.

The right-hand side of vorticity (3.4), (3.7) for VE fluids in addition to the gravity term includes ${\boldsymbol {\nabla }}\times {({\boldsymbol {\psi }}/{\rho _d})}$. This term incorporates the TS waves. In our previous work (Dharodi et al. Reference Dharodi, Tiwari and Das2014, Reference Dharodi, Das, Patel and Kaw2016), in absence of a gravity term, we observed such waves emerging from the rotating vorticity lobes into the medium.

3.6. Rising bubble

In order to simulate the case of a rising bubble, the numerical simulations have been carried out for ${a_c}=2.0$, ${\bar {\rho }_{d1}}=-0.5$, ${x_{c}}=0$ and ${y_{c}}=-4{\rm \pi}$ in (3.11). The homogeneous background density ${\rho _{d0}}=5$. Figure 8 displays the time evolution of the bubble density profile of an inviscid fluid. Initially, at time $t=0$, the bubble is axisymmetric (see figure 6a) and as the system evolves the bubble simply rises without any significant change in shape for a short period of time. At a later stage, the initially circular profile assumes a crescent shape, as evident from the first panel. Further, as time progresses, a mushroom-like structure that is characteristic of RT instability begins to appear along with rolling and mixing at the edge of the bubble. This mushroom-like structure then breaks into two distinct elliptical blobs as evident from the figure. These blobs propagate forward as a single entity leaving behind the wake-like structure in background fluid.

Figure 8. Time evolution of bubble density (see figure 11 for the respective vorticity) profile for inviscid HD fluid.

Next, with the understanding of bubble evolution for inviscid fluids, the time evolution of bubble for VE fluids has been performed. To envisage the effect of coupling strength, we compare one VE fluid (figure 9a–d) having ${\eta }=2.5, {\tau _m}=20$, i.e. $\eta /\tau _m=0.125$ with another VE fluid (figure 9e–h) having higher coupling strength, i.e. ${\eta }=5$, ${\tau _m}=10$, $\eta /\tau _m=0.5$. We observe a decrease in the vertical upward motion of the bubble and larger separation between elliptical blobs in the horizontal direction earlier in figure 9(e–h) than in figure 9(a–d). At the same time, the rolling rate of blobs increases causing the expansion of blobs and the disappearance of the dragging tail. Thus, we conclude that similar to the RT instability, the growth of the bubble (vertical upward motion) gets reduced with increasing coupling strength. The cause of this reducing effect can be easily understood through the respective vorticity evolution given in figure 10. In figure 10, the vorticity plots in panels (a–d) and (e–h) correspond to the density profiles shown in figures 9(a–d) and 9(e–h), respectively.

Figure 9. Evolution of bubble density profile (see figure 10 for the respective vorticity evolution) in time for two VE fluids: (a–d) $\eta =2.5$, $\tau _m=20$ (see supplemental material, https://doi.org/10.1017/S0022377821000349) and (e–h) $\eta =5$, $\tau _m=10$. The growth of the bubble (vertical upward motion) gets reduced with increasing coupling strength.

Figure 10. Time evolution of bubble vorticity profile (see figure 9 for the respective density evolution) for two VE fluids: (a–d) $\eta =2.5$, $\tau _m=20$ (see supplemental material) and (e–h) $\eta =5$, $\tau _m=10$. There is emission of a TS wave surrounding the vorticity lobes whose phase velocity is proportional to the coupling strength. It is clear that the increase of phase velocity of TS waves with coupling strength suppresses the BD instability (i.e. reduction in vertical upward motion of lobes).

It is evident from figure 10 that there is emission of TS waves surrounding the vorticity lobes for VE fluids. The emission of waves from either lobe pushes the other lobe in the direction perpendicular to the direction of propagation of the entire structure and as a result the lobes get well separated with time and proportional to the coupling strength (a detailed discussion on the evolution of such dipolar vorticity structures in VE fluids is given by Dharodi et al. (Reference Dharodi, Das, Patel and Kaw2016)). Besides this lobe separation, the emission of TS waves notably reduced the strength of the dipole thereby delaying the dipole propagation and this reduction in dipole propagation is also proportional to the coupling strength. The relative observations of figure 10(a–d) and figure 10(e–h) clearly reflect it through the TS enclosure which is larger for the $\eta =5$, $\tau _m=10$ than $\eta =2.5$, $\tau _m=20$. This reduction in vertical propagation direction and enhancement in horizontal separation between lobes/blobs with increasing coupling strength means the suppression of the buoyant nature of a bubble is proportional to the coupling strength of the medium. There is no such TS waves surrounding the vorticity for inviscid HD fluid (figure 11).

Figure 11. Evolution of bubble vorticity profile in time for inviscid HD fluid (see figure 8 for the respective density).

To substantiate our observations made so far, four snapshots of density and vorticity evolution of bubble at the same time step $({\approx }20)$ for four different combinations of $\eta$ and $\tau _m$ are shown in figures 12 and figure 13, respectively. From figure 12, it is clear that the BD instability gets weaker with increasing $\eta$ (in the vertical direction with a fixed value of $\tau _m$), while $\tau _m$ speed up the instability (along the horizontal direction with a fixed value of $\eta$). Snapshots also show that the instability evolution is similar for the same ratio ${\eta /\tau _m}$. The role of TS waves can be visualized from the vorticity snapshots given in figure 13.

Figure 12. Plots of $\eta$ versus $\tau _m$. The snapshots of bubble density profiles at a particular time step (see figure 13 for the respective vorticity). The BD instability gets weaker with increasing $\eta$ (in the vertical direction with a fixed value of $\tau _m$) and $\tau _m$ speeds up the instability (along the horizontal direction with a fixed value of $\eta$).

Figure 13. Plots of $\eta$ vs $\tau _m$. The snapshots of bubble vorticity profiles at a particular time step (see figure 12 for the respective density). The BD instability gets weaker with increasing viscosity (in the vertical direction with a fixed value of $\tau _m$) and $\tau _m$ speeds up the instability (along the horizontal direction with a fixed value of $\eta$).

3.7. Falling droplet

For simulating the falling droplet, i.e. case B (figure 6b), the values of parameters in (3.11) for the droplet density profile are ${a_c}=2.0$, ${\bar {\rho }_{d1}}=0.5$, ${x_{c}}=0$, ${y_{c}}=4{\rm \pi}$ and ${\rho _{d0}}=5$. The total density is $\rho _d=\rho _{d0}+\rho _{d1}$. Figure 14 shows the dynamics of a falling droplet in tranquil inviscid HD fluid. As the simulation begins, under the influence of gravity, the suspended drop starts to fall. As time passes, the drop breaks up forming first a semicircular structure then two blobs similar to the case of bubble. Figure 15 reveals the different stages of droplet for VE fluids. As done for the rising bubble, we shall now compare figure 15(a–d) having ${\eta }=2.5, {\tau _m}=20$, i.e. $\eta /\tau _m=0.35$ with figure 15(e–h) having ${\eta }=5, {\tau _m}=10$, i.e. $\eta /\tau _m=0.5$. For a higher coupling strength shown in figure 15(e–h), the downward motion is slower. The transverse dimension is even larger, and the dragging tail does not seem to appear at all. Again, similar to the case of bubbles, the vorticity plots are provided for observing a clear role of TS waves, i.e. pushing the two lobes apart, as seen in figure 16. For the inviscid HD fluid, the corresponding vorticity subplots (figure 17) to the mentioned density profile (figure 14), no such wave are observed which could constraint the vertical downward motion. Thus, we conclude that the increasing ratio $\eta /\tau _{m}$ suppresses the buoyant nature of a droplet.

Figure 14. Time evolution of droplet density (see figure 17 for the respective vorticity evolution) profile for inviscid HD fluid.

Figure 15. Evolution of droplet density profile (see figure 16 for the respective vorticity evolution) in time for two VE fluids: (a–d) $\eta =2.5$, $\tau _m=20$ (see supplemental material) and (e–h) $\eta =5$, $\tau _m=10$. The growth of the droplet (vertical downward motion) gets reduced with increasing coupling strength.

Figure 16. Time evolution of droplet vorticity profile (see figure 15 for the respective density evolution) for two VE fluids: (a–d) $\eta =2.5$, $\tau _m=20$ (see supplemental material) and (e–h) $\eta =5$, $\tau _m=10$. There is emission of a TS wave surrounding the vorticity lobes whose phase velocity is proportional to the coupling strength. It is clear that the increase of phase velocity of TS waves with coupling strength suppresses the BD instability (i.e. reduction in vertical downward motion of lobes).

Figure 17. Time evolution of droplet vorticity (see figure 14 for the respective density evolution) profile for inviscid HD fluid.

In the same way as the bubble evolution, figures 18 and figure 19 show the snapshots of droplet density and vorticity evolution at a time step=$({\approx }20)$, respectively. Here, also, from figure 18, it is clear that the BD instability gets weaker with increasing $\eta$, while $\tau _m$ speed up the instability. Snapshots, also, show that the instability evolution is similar for the same ratio ${\eta /\tau _m}$. The role of TS waves can be visualized from the vorticity snapshots given in figure 19.

Figure 18. Plots of $\eta$ vs $\tau _m$. The snapshots of droplet density profiles at a particular time step (see figure 19 for the respective vorticity). The BD instability gets weaker with increasing viscosity (in the vertical direction with a fixed value of $\tau _m$) and $\tau _m$ speeds up the instability (along the horizontal direction with a fixed value of $\eta$).

Figure 19. Plots of $\eta$ vs $\tau _m$. The snapshots of droplet vorticity profiles at a particular time step (see figure 18 for the respective density). The BD instability gets weaker with increasing viscosity (in the vertical direction with a fixed value of $\tau _m$) and $\tau _m$ speeds up the instability (along the horizontal direction with a fixed value of $\eta$).