2.1 General description

We consider a spherical droplet of diameter $D=2R$ , of a liquid of density $\unicode[STIX]{x1D70C}_{l}$ and viscosity $\unicode[STIX]{x1D707}_{l}$ falling from a small height $h$ between the droplet and the substrate, with an initial velocity $U_{0}$ normal to the substrate, surrounded by a gas of density $\unicode[STIX]{x1D70C}_{g}$ , viscosity $\unicode[STIX]{x1D707}_{g}$ , the liquid–gas surface tension being $\unicode[STIX]{x1D6FE}$ (see figure 2). The substrate is dry and flat and a typical partially wettable surface is considered with a contact angle $\unicode[STIX]{x1D703}=90^{\circ }$ for the sake of simplicity. We have also tested other contact angles (from $\unicode[STIX]{x1D703}=0$ to $\unicode[STIX]{x1D703}=180^{\circ }$ ) finding no qualitative changes in the mechanisms described below, so that we only show here the results for $\unicode[STIX]{x1D703}=90^{\circ }$ . It also indicates that the details of the moving contact line does not play a significant role in the impact dynamics, because of the rapid spreading, although it could be crucial in the receding dynamics of the drop as observed experimentally (Bartolo, Josserand & Bonn Reference Bartolo, Josserand and Bonn2005). Notice however, that the full account of the moving contact-line dynamics can be implemented in numerical simulations using specific boundary conditions (Afkhami, Zaleski & Bussmann Reference Afkhami, Zaleski and Bussmann2009; Legendre & Maglio Reference Legendre and Maglio2015; Mahadi, Afkhami & Kondic Reference Mahadi, Afkhami and Kondic2015; Mahady, Afkhami & Kondic Reference Mahady, Afkhami and Kondic2016; Afkhami et al. Reference Afkhami, Buongiorno, Guion, Popinet, Scardovelli and Zaleski2017). Experiments are typically achieved by releasing a pendant liquid droplet from a height $H\gg D$ , that then falls under gravity, so that $U_{0}\sim \sqrt{gH}$ . If the Froude number ( $Fr=U_{0}^{2}/(gD)\sim H/D$ ) is large, the liquid inertia dominates the gravitational force so that gravity can be neglected during impact itself. We assume a spherical droplet although its shape can sometimes be affected both by the initial oscillations triggered by the dynamics of detachment (Wang, Kuan & Tsai Reference Wang, Kuan and Tsai2013; Thoraval, Li & Thoroddsen Reference Thoraval, Li and Thoroddsen2016) and the gas drag (Thoroddsen et al. Reference Thoroddsen, Etoh, Takehara, Ootsuka and Hatsuki2005). Indeed, we will focus here on the short time dynamics that will be mostly determined by the radius of curvature of the bottom of the drop that we suppose to be close to $R$ . In general, drop impact is characterized by two dimensionless numbers based on the liquid drop parameters: the Reynolds number that compares inertial effects with viscous effects

and the Weber number which compares inertial effects with capillary effects

Here, since we neglect gravity, compressibility and non-continuum effects, only two additional dimensionless numbers are present in the problem: the density and viscosity ratio $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D70C}_{g}/\unicode[STIX]{x1D70C}_{l}$ and $m=\unicode[STIX]{x1D707}_{g}/\unicode[STIX]{x1D707}_{l}$ respectively. Note that the incompressible lubrication of the gas layer is controlled by the Stokes number alone:

which is nothing more than the ratio between the viscosity ratio and the Reynolds number. Compressibility is often neglected for droplet impact problems as long as the impact velocity is much below the speed of sound in both fluids (Lesser & Field Reference Lesser and Field1983). In fact, if this assumption is usually correct for most of the flow and in the liquid domain, it has been shown that it is not the case in the small region beneath the drop just before the impact if the gas pressure is small enough (Mandre et al. Reference Mandre, Mani and Brenner2009). There, because of the high pressure coming from the inertia of the drop cushioning the thin gas film, the lubrication gas pressure can overcome the ambient pressure so that compressible effects have to be taken into account. The potential influence of compressibility on the splashing dynamics has been studied in great detail in the framework of the simplified model solving the lubrication equation for the gas and an inviscid dynamics for the liquid drop approximated by a parabolic shape (Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010; Mandre & Brenner Reference Mandre and Brenner2012). Although it shows that changing the surrounding gas pressure, as done in the pioneering experiments of Xu et al. (Reference Xu, Zhang and Nagel2005), can lead to different compressible regimes, no clear relation between this dynamics and the suppression of the splash was deduced beside qualitative arguments (Mandre & Brenner Reference Mandre and Brenner2012). This is why in the present paper we propose to investigate the dynamics omitting the compressibility of the gas in order to isolate the different effects leading to the splash in the dynamics. In particular, our approach will help to identify a purely incompressible inertial contribution favouring the splashing at high gas pressure (in our case it corresponds then to high gas density), explaining thus how a lower gas density could eventually suppress the splash.

Figure 2.

Calculation domain. The domain is constituted of a square. A semi-circular droplet of diameter $D$ , density $\unicode[STIX]{x1D70C}_{l}$ and viscosity $\unicode[STIX]{x1D707}_{l}$ falls from a height $h$ with an initial velocity $U_{0}$ . The surrounding gas has a density $\unicode[STIX]{x1D70C}_{g}$ , a viscosity $\unicode[STIX]{x1D707}_{g}$ and a surface tension $\unicode[STIX]{x1D6FE}$ . The bottom boundary ( $z=0$ ) is considered as a solid substrate with no-slip condition and a contact angle $\unicode[STIX]{x1D703}=90^{\circ }$ is defined.

Within this framework, the system obeys the Navier–Stokes equations for incompressible flow with surface tension and sharp interfaces that read:

where $p$ and $\boldsymbol{u}$ are the pressure and velocity fields respectively. $\mathbb{D}=1/2(\text{}^{t}\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u})$ is the rate-of-strain tensor ( $\mathbb{D}_{ij}\equiv (\unicode[STIX]{x2202}_{i}u_{j}+\unicode[STIX]{x2202}_{j}u_{i})/2$ ). The distribution function $\unicode[STIX]{x1D6FF}_{s}$ is a delta function which insures the jump conditions at the interface, $\unicode[STIX]{x1D705}$ the curvature and $\boldsymbol{n}$ the normal vector to the interface. In this formulation the density $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)$ and viscosity $\unicode[STIX]{x1D707}(\boldsymbol{x},t)$ are fields that take constant values ( $\unicode[STIX]{x1D70C}_{l}$ and $\unicode[STIX]{x1D70C}_{g}$ for the density, $\unicode[STIX]{x1D707}_{l}$ and $\unicode[STIX]{x1D707}_{g}$ for the viscosity). Finally, since we are interested in the short-time dynamics after impact, we can make the assumption that the dynamics is axisymmetric, although the influence of three-dimensional effects would need to be investigated in the future.

2.2 Numerical method

The set of incompressible two-fluid Navier–Stokes equations (2.4) and (2.5) is solved in dimensionless form using the free code Gerris (Popinet Reference Popinet2003, Reference Popinet2006, Reference Popinet2009) in a two-dimensional axisymmetric geometry (using cylindrical space variable $r$ and $z$ ). Rescaling the density, space, time, velocity and pressure by $\unicode[STIX]{x1D70C}_{l}$ , $D$ , $D/U_{0}$ , $U_{0}$ and $\unicode[STIX]{x1D70C}_{l}U_{0}^{2}$ respectively, the Navier–Stokes equations become (using the $^{\ast }$ for dimensionless variables):

involving the Reynolds and Weber numbers defined above. The density $\unicode[STIX]{x1D70C}^{\ast }$ and viscosity $\unicode[STIX]{x1D707}^{\ast }$ fields now take values $1$ in the liquid phase and $\unicode[STIX]{x1D6FC}$ and $m$ in the gas phase. The dimensionless strain tensor now reads $\mathbb{D}_{ij}^{\ast }\equiv (\unicode[STIX]{x2202}_{i}^{\ast }u_{j}^{\ast }+\unicode[STIX]{x2202}_{j}^{\ast }u_{i}^{\ast })/2$ . For the sake of simplicity, we will drop the $^{\ast }$ in what follows. Gerris solves this set of equations using a finite-volume method on an adaptive tree-structured (quadtree/octree) mesh. Momentum and pressure are computed by a projection method and the velocity advection is calculated by a second-order Bell–Collela–Glaz scheme. The interface is tracked by a volume-of-fluid/piecewise linear interface calculation (VOF/PLIC) method (Li Reference Li1995) with a mixed Young’s centred scheme (Tryggvason, Scardovelli & Zaleski Reference Tryggvason, Scardovelli and Zaleski2011) to determine the normal vector and a Lagrangian-explicit scheme for VOF advection. Curvature is computed by the height-function method and surface tension is estimated from curvature by a balanced continuous-surface-force method (Popinet Reference Popinet2009). The density and viscosity fields are computed based on the VOF fraction $C$ using an arithmetic mean. The adaptive mesh refinement is controlled using the following criteria, as illustrated in figure 3:

1. (i) the gradient of the VOF function $C$ tracking the interface,

2. (ii) the variation of the vorticity,

3. (iii) the curvature in the meshes near the interface,

4. (iv) the solid surface.

The refinement consists in dividing a square parent cell into four (eight in three dimensions) square children cells of size half that of the parent one. Therefore, the level of refinement $n$ corresponds to a cell size $2^{-n}$ smaller than the computational unit size. Finally, a no-slip boundary condition is imposed on the solid substrate boundary ( $z=0$ ), while free slip is imposed on the other two boundaries (the last one, corresponding to $r=0$ comes from the axial symmetry). Because of the Navier–Stokes solver, this no-slip boundary condition imposes in fact an effective slip length of the order of half of the grid size at the substrate and a convergence study needs to be performed to investigate the influence of such a condition on the results.

Figure 3.

Snapshots of the adaptive mesh refinement at two different instants of the impact: (a) just before impact where most of the refinement is concentrated on the interface and solid surface, and (b) a zoom after the impact when a layer of liquid is skating above the solid surface, the red curve representing the interface.

In this paper, since we want to characterize the influence of the gas properties (density and viscosity) on the dynamics, we have kept all the other quantities constant, specifically $Re=1000$ and $We=370$ . In other words we are performing a parametric study by varying only the ratios $\unicode[STIX]{x1D6FC}$ and $m$ (notice that the Stokes number also varies since the gas viscosity does). These values of the Reynolds and Weber numbers have to be compared with those of the original experiment of Xu et al. (Reference Xu, Zhang and Nagel2005), $Re\sim 10\,000$ and $We\sim 1500$ that are significantly higher. In fact, numerical simulations at such high Reynolds and Weber numbers are possible but may rapidly suffer from a lack of spatial resolution because of the thin liquid layer formed at impact. The influence of the gas density is however present for the values of the parameters chosen here so that a qualitative parametric study of the mechanisms involved is possible. Furthermore, in the actual configuration, a typical simulation run takes about 12 hours on one core of a powerful machine with a refinement level $n=11$ while it would take almost $5$ days for $n=12$ , making the parametric study done here difficult. Finally, our numerics can be understood as the impact of a drop of smaller diameter than the original experiment, with a higher liquid viscosity, smaller velocity and smaller surface tension. For instance, considering the velocity and the drop diameter both divided by two compared to the experiment, one would obtain our numerical value by multiplying the liquid viscosity by 2.5 and dividing the surface tension by 2, something that can easily be obtained experimentally. The initial distance $h$ between the droplet and the substrate has been set such that it has no influence on the impact dynamics, leading to $h=0.1$ in dimensionless unit. A typical impact obtained by numerical simulations is shown in figure 1 for $\unicode[STIX]{x1D6FC}=0.003$ and $m=0.15$ . The origin of time is taken as the virtual time at which a perfectly spherical droplet falling at constant speed would contact the solid providing the droplet does not deform. In other words, the numerical simulation starts at dimensionless time $t=-h$ . As the droplet approaches the solid substrate, the gas cushioning delays the contact between the droplet and the substrate that happens slightly later than for $t=0$ and not at the droplet bottom ( $r=0$ ) but rather at a positive radius: a small bubble of gas is subsequently entrapped by the dynamics between the expanding liquid drop and the solid substrate. This bubble cannot really be distinguished in figure 1 but its formation can be observed on the refined mesh of figure 3(b) which is taken before the bubble entrapment, and will be characterized below. The drop then spreads on the substrate and a thin jet emerges at some angle, forming the splash. The splash deforms with time and a characteristic S-shape is created (Josserand et al. Reference Josserand, Ray and Zaleski2016).

Numerical convergence for drop impact simulation is a hard challenge because of the small structures created by the impact (Josserand et al. Reference Josserand, Ray and Zaleski2016) but also due to the moving contact-line dynamics that has to be accounted for in the present configuration. The numerical convergence of our results has been tested by measuring the minimal film thickness, which is one of the most sensitive quantities among those that depend on the mesh size and on the boundary condition on the substrate. Figure 4 shows this quantity as a function of time for increasing mesh refinement demonstrating a good convergence for refinement levels above $n=10$ . For the rest of this paper, all the results will be computed using the same maximum refinement level $n=11$ . However, the convergence for the dynamics of the moving contact line cannot be obtained in our numerics because of the no-slip boundary conditions that leads to a singular behaviour in the continuum limit (Dussan Reference Dussan1979). Physically, the flow is regularized by small-scale effects, that can be accounted for through an effective boundary slip condition at the solid surface, with a slip length of the order of a few tens of nanometres (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). To obtain numerical convergence for this dynamics, complex boundary conditions at the moving contact line involving the mesh size have to be introduced (Afkhami et al. Reference Afkhami, Zaleski and Bussmann2009; Legendre & Maglio Reference Legendre and Maglio2015; Afkhami et al. Reference Afkhami, Buongiorno, Guion, Popinet, Scardovelli and Zaleski2017). Keeping a fixed maximum refinement level avoids the question of the numerical convergence of the moving contact-line dynamics, since it corresponds then to a fixed slip length boundary condition at the solid substrate. For instance, for a millimetre drop diameter, the slip length in our simulations is approximately 250 nanometres, clearly larger than the physical one, but with no expected qualitative change in the dynamics.

Figure 4.

Evolution with time of $h_{min}(t)$ , the minimal air thickness in the numerics, for different maximum refinement levels $n$ (numbers indicated bottom left of the figure) for $m=0.1$ and $\unicode[STIX]{x1D6FC}=0.001$ . The grid size is then proportional to $1/2^{n}$ near the substrate. Although small deviations can be observed between the different levels, we observe a rapid convergence of the results. The straight line indicates the free fall of the drop and the inset shows a zoom of the figure near the time where the bubble is entrapped.

2.3 Scaling analysis

Prior to performing a parametric study of the dynamics, it is useful to define and to estimate some quantities related to the impact. In particular, a geometrical argument based on the intersection of the impacting drop with the substrate shows that, at short times, the vertical length of the intersection $h(t)$ varies like $t$ while the horizontal one follows $r_{c}(t)\sim \sqrt{h}\sim \sqrt{t}$ . This simple argument has been used since Wagner (Reference Wagner1932), implying a self-similar dynamics of the impact at short times (Josserand & Zaleski Reference Josserand and Zaleski2003; Smith et al. Reference Smith, Li and Wu2003). More precisely, a quantitative estimate can even be done for the foot of the jet following $\sqrt{3t/2}$ (Riboux & Gordillo Reference Riboux and Gordillo2015; Philippi, Lagrée & Antkowiak Reference Philippi, Lagrée and Antkowiak2016) (notice that the $3/2$ factor, in contrast with the 3 in the literature, comes from the rescaling by $D$ used here instead of $R$ ) but remarkably, the general square-root-of-time law (with possibly a different prefactor) describes correctly all the horizontal scales at small times (Josserand & Zaleski Reference Josserand and Zaleski2003). Figure 5 shows the spreading radius $r_{w}(t)$ , defined as the position where the interface height is minimum, as a function of time for $\unicode[STIX]{x1D6FC}=0.001$ and $m=0.1$ as an illustration of this evolution. In fact, it has been shown that the spreading radius varies only slightly with the fluid viscosity and surface tension (Josserand & Zaleski Reference Josserand and Zaleski2003), so that the curves obtained with other parameters will show the same overall aspect. The curve is well fitted by the predicted law:

where $t_{d}$ is the time when the dimple starts to form beneath the drop. Remarkably, the fit overestimates the spreading radius at very short time, probably because of a small capillary correction when the dimple forms. Rapidly however, the predicted law gives a very good approximation of the spreading radius.

Figure 5.

Radius $r_{w}$ , defined as the radius of minimal height as a function of time $t-t_{d}$ , where $t_{d}=0.102$ is the time at which the dimple beneath the drop appears, for $\unicode[STIX]{x1D6FC}=0.001$ and $m=0.1$ . The inset shows the relation between $r_{w}$ and $\sqrt{t-t_{d}}$ . Symbols denote the temporal evolution of the $r$ -coordinate of the $h_{min}$ point and the dashed curve represents the law $r_{w}=\sqrt{3(t-t_{d})/2}$ .

Small discontinuities can be observed above $t-t_{d}\sim 0.03$ corresponding to $r_{w}\sim 0.2$ . Above these values, the central bubble has been entrapped and the dynamics of the numerical contact line induces a sequence of contact zone separated by bubbles as illustrated in figure 6. There, $r_{w}(t)$ is defined as the largest value of the liquid in contact with solid, explaining the discontinuities in the curve 5. These bubbles are formed by an instability of the moving contact line but note that due to the cylindrical symmetry of our calculations these bubbles are in fact tori that would be unstable in three dimensions. These bubbles are similar to those observed in numerical simulations of drop impacts on a thin liquid film (Josserand et al. Reference Josserand, Ray and Zaleski2016) and reminiscent of experimental observations by Li, Vakarelski & Thoroddsen (Reference Li, Vakarelski and Thoroddsen2015b ). How these bubbles would fare under refinement of the grid has not been investigated but their size would likely vary.

Figure 6.

The expanding lamella after bubble entrapment ( $t=0.059$ ) for the same parameters as in figure 5. The liquid–gas interface is drawn while the solid surface is the dashed line plotted at $z=0$ . Small bubbles are entrapped within the wetted area due to the rapid expansion of the lamella, as displayed in the inset that shows a zoom of the wetted area.

This well-known bubble entrapment mechanism has been observed experimentally (Thoroddsen, Etoh & Takehara Reference Thoroddsen, Etoh and Takehara2003; Thoroddsen et al. Reference Thoroddsen, Etoh, Takehara, Ootsuka and Hatsuki2005) and numerically (Mehdi-Nejad et al. Reference Mehdi-Nejad, Mostaghimi and Chandra2003; Smith et al. Reference Smith, Li and Wu2003; Korobkin et al. Reference Korobkin, Ellis and Smith2008; Mandre et al. Reference Mandre, Mani and Brenner2009; Duchemin & Josserand Reference Duchemin and Josserand2011) in many studies and contexts. It is usually explained using the lubrication equation framework in the gas (Smith et al. Reference Smith, Li and Wu2003; Korobkin et al. Reference Korobkin, Ellis and Smith2008; Mandre et al. Reference Mandre, Mani and Brenner2009; Mani et al. Reference Mani, Mandre and Brenner2010) and we will rephrase now the scaling arguments for consistency. Indeed, when the drop approaches the solid, the gas layer beneath the drop becomes very thin so that the gas inertia can be neglected in the Navier–Stokes equations leading to the lubrication equation for the gas layer thickness $h(r,t)$ , involving the pressure $p_{g}(r,t)$ in the gas, written here in dimensionless form:

Using the scalings induced by the geometrical argument above and considering that $\unicode[STIX]{x2202}_{t}h\sim 1$ , we obtain the following scaling for the lubrication pressure in the gas:

On the other hand the pressure in the liquid $p_{l}$ needed for the deformation of the drop at impact can be estimated using the same self-similar argument by balancing the pressure force $p_{l}\times \unicode[STIX]{x03C0}r_{c}^{2}$ with the momentum variation in the half-sphere of radius $r_{c}$ (in dimensionless form, the vertical velocity is $U^{\ast }=1$ , and the liquid density is $\unicode[STIX]{x1D70C}^{\ast }=1$ , kept in the first formula below to avoid confusion):

yielding (taking now $\unicode[STIX]{x1D70C}^{\ast }=U^{\ast }=1$ and using the fact that $r_{c}(t)=\sqrt{t}$ so that $\text{d}r_{c}/\text{d}t=1/(2r_{c})$ )

which leads to the scaling for the pressure in the liquid near the impact zone:

The entrapment of the bubble can be interpreted within this simplified scaling analysis: when $t$ is small enough, the lubrication pressure $p_{g}$ is much bigger than the inertial one $p_{l}$ needed to deviate the impacting drop so that the drop skates on a thin gas layer. On the other hand, when $t$ increases, the lubrication pressure decreases much more rapidly than the inertial one so that at a critical time, the gas layer cannot any more deviate the liquid that then enters in contact with the solid. This threshold, measured in term of entrapment time $t_{e}$ or height $h_{e}=t_{e}$ , radius $r_{e}=\sqrt{t_{e}}$ and volume $V_{e}\propto r_{e}^{2}h_{e}$ can be estimated by balancing the two pressures $p_{l}\sim p_{g}$ leading to:

Measuring these quantities experimentally remains a difficult challenge because of the small values of the Stokes numbers but this has been done eventually using complex optical techniques (Driscoll & Nagel Reference Driscoll and Nagel2011; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Kolinski et al. Reference Kolinski, Rubinstein, Mandre, Brenner, Weitz and Mahadevan2012; Klaseboer, Manica & Chan Reference Klaseboer, Manica and Chan2014; Li, Vakarelski & Thoroddsen Reference Li, Vakarelski and Thoroddsen2015a ). The results are somehow puzzling since the theoretical predictions have not always been observed or only in a narrow range of Stokes numbers. We have investigated numerically the radius of the entrapped bubble by measuring the radius $r_{e}$ of the first contact between the droplet and the solid and the height $h_{e}$ of the gas film at that time taken at $r=0$ , for different Stokes numbers, by varying the gas viscosity. As shown in figure 7(a), the scaling laws work well for a large range of Stokes numbers, proving that the lubrication is clearly in action at first order. However they fail for small or very large Stokes numbers: the finite size of the drop affects the scaling for high Stokes numbers while one may expect inertia effects to dominate lubrication theory at small Stokes numbers (Wilson & Duffy Reference Wilson and Duffy1998; Luchini & Charru Reference Luchini and Charru2010), as already observed numerically for the impact on a liquid surface (Jian et al. Reference Jian, Josserand, Ray, Duchemin, Popinet and Zaleski2015; Josserand et al. Reference Josserand, Ray and Zaleski2016). This variation of the entrapped bubble size with the gas density at small Stokes numbers is a first manifestation of an effect related to the gas inertia in the impact. This is shown in figure 7(b), where the volume is plotted as the density ratio varies for a fixed Stokes number ( $St=10^{-5}$ ), exhibiting a decrease of the bubble entrapped with decreasing gas density. There, while a monotonic increasing trend can be seen, a high scattering of the data is also observed. This scattering comes from the small gas density that enhances numerical instability. In particular, large velocity gradients form in this region where the maximum refinement level is reached, so that these gradients are probably not sufficiently well resolved. Future studies should use higher refinement to clarify this point. When comparing the lubrication pressure with the Bernoulli pressure in this configuration, it has been shown that the threshold between inertia and lubrication occurs for $\unicode[STIX]{x1D6FC}\sim St^{2/3}$ (Josserand et al. Reference Josserand, Ray and Zaleski2016), in good qualitative agreement with the plateau observed at small density ratio in figure 7(b). Even if it is tempting to associate this volume decrease with a smaller angle of jet formation and thus with a weaker splash, the full dynamics needs to be investigated further to elucidate the splashing mechanisms.

Figure 7.

(a) Thickness of the gas layer $h_{e}$ , radius of the first contact point $r_{e}$ and volume of the central entrapped bubble $V_{e}$ versus the Stokes number $St$ (logarithmic scale) for density ratio $\unicode[STIX]{x1D6FC}=0.003$ (this implies  $m=Re\cdot St=1000\,St$ ), symbols represent simulation results. Dashed lines represent the scaling laws $h_{e}=2.6\,St^{2/3}$ , $r_{e}=3.3\,St^{1/3}$ and $V_{e}=28.314St^{4/3}$ respectively. (b) $V_{e}$ versus density ratio $\unicode[STIX]{x1D6FC}$ for $St=10^{-5}$ .