23.1 Introduction

We will now turn to an effect that is commonly observed whenever a thin thread of liquid is running, e.g., from the tap. The thread of the liquid will remain stable for some time before eventually decomposing into individual droplets. This effect is referred to as Plateau-Rayleigh^1,2instability. A fluid column experiences significant forces originating from surface tension. These can be minimized when minimizing the surface area. Therefore, a liquid column will have the tendency to decompose into individual droplets. In this section, we will study the fluid mechanics of this phenomenon.

23.2 Stability Considerations

In general, the instability of a thin cylinder of a liquid is due to that fact it usually has a higher surface compared to a sphere. We will discuss a thin fluid cylinder with radius r_cyl and length l_cyl. Its surface-to-volume ratio is given by

${\frac{A}{V}|}_{\text{cyl}}=\frac{2\pi r_{\text{cyl}}^2+2\pi r_{\text{cyl}}{l}_{\text{cyl}}}{\pi r_{\text{cyl}}^2{l}_{\text{cyl}}}2\left(\frac{1}{l_{\text{cyl}}}+\frac{1}{r_{\text{cyl}}}\right)$

For a sphere of diameter r_sphere the surface-to-volume ratio is

${\frac{A}{V}|}_{\text{sph}}=\frac{4\pi r_{\text{sphere}}^2}{\frac{4}{3}\pi r_{\text{sphere}}^2}=\frac{3}{r_{\text{sphere}}}$

In order to determine if it is energetically more favorable to split the cylinder into individual droplets, the ratio of the two surface-to-volume ratios has to be found. Using Eq. 23.1 and Eq. 23.2 we find

$\begin{array}{ll}\frac{{\frac{A}{V}|}_{\text{cyl}}}{{\frac{A}{V}|}_{\text{sph}}}\hfill & =\frac{2\left(\frac{1}{l_{\text{cyl}}}+\frac{1}{r_{\text{cyl}}}\right)}{\frac{3}{r_{\text{sphere}}}}=\frac{2}{3}\left(\frac{r_{\text{sphere}}}{l_{\text{cyl}}}+\frac{r_{\text{sphere}}}{r_{\text{cyl}}}\right)\hfill \\ \hfill & \approx \frac{2}{3}\frac{r_{\text{sphere}}}{l_{\text{cyl}}}\hfill \end{array}$

where we used the last approximation because r_sphere ˂˂ l_cyl. From Eq. 23.3 we can see that for $r_{\text{sphere}}>\frac{2}{3}{r}_{\text{cyl}}$ it is energetically more favorable for the fluid cylinder to decompose into individual droplets. Whether or not this will happen will depend on the inertia of the fluid, i.e., the forces that hold the cylinder compact, the surface tension, and the gravitational forces that may accelerate and thus stretch the cylinder. These fluid mechanical considerations will be discussed in the following section at the falling fluid jet.

23.3 Fluid Jets

The simplest form of a liquid jet is the falling fluid jet, which essentially is a falling cylinder of a given fluid. This fluid exits an opening at a given initial velocity and thus has a certain inertia that keeps it compact. The cylinder is accelerated by gravity and experiences surface tension forces. Fig. 23.1 depicts the scenario of this falling fluid jet that we know, e.g., from tap water. In the following we will derive the balance of forces at the fluid jet. The two parameters we want to derive are the average velocity v_z (z) and the radius r (z) at the location z. Please note that the independent variable z is chosen along the path of the falling jet, i.e., in the direction of the gravitational force. Besides gravity, there are other forces that may be used to drive a fluid jet, e.g., viscous forces used, e.g., in hydrodynamic focusing or droplet microfluidics.

Before we start with the derivation of the fluid mechanics, we will shortly reflect on what we expect to happen in a qualitative manner. The liquid jet is created by pressing a liquid through an opening, i.e., by means of a nozzle structure. The liquid will exit with a given inertia. This inertia will tend to keep the liquid jet compact, i.e., in form of a cylinder. While falling, the liquid will be accelerated, which, as discussed, will lengthen it. This in consequence should result in a decrease of the jet’s radius. This will induce curvatures at the jet’s surface. These curvatures will cause surface tension that, as discussed, will try to minimize the surface of the jet, eventually leading to its decomposition into droplets.

We are not interested in the velocity profile across the radius v_z (r, z), but will regard the velocity in the z-direction to be only a function of the independent variable z. From Bernoulli’s equation (see Eq. 14.11) we know that

$\frac{\rho }{2}{\upsilon }_{z,0}^2+\rho g{z}_0+p_1=\frac{\rho }{2}{\upsilon }_z^2\left(z\right)+\rho gz\left(z\right)+p\left(z\right)$

From Eq. 20.11 we know that the respective pressures can be expressed as functions of the curvature of the surface of the jet

$\begin{array}{ll}{p}_1-p_1\hfill & =\gamma \left(\frac{1}{r_0}+\frac{1}{r_{\infty }}\right)\hfill \\ \hfill & =\frac{\gamma }{r_0}\hfill \end{array}$

$\begin{array}{ll}p\left(z\right)-p_1\hfill & =\gamma \left(\frac{1}{r\left(z\right)}+\frac{1}{r_{\infty }}\right)\hfill \\ \hfill & =\frac{\gamma }{r\left(z\right)}\hfill \end{array}$

where p₁ is the ambient pressure. The term $\frac{1}{r_{\infty }}$ can be set to 0 as the jet would (in the unperturbed state) be a straight cylinder, and therefore the curvature of its surface along the z-axis would be 0. Using Eq. 23.5 and Eq. 23.6 we can rewrite Eq. 23.4 as

$\begin{array}{l}\frac{\rho }{2}{\upsilon }_{z,0}^2+\rho g{z}_0+p_1-p\left(z\right)=\frac{\rho }{2}{\upsilon }_z^2\left(z\right)+\rho gz\left(z\right)\\ \frac{\rho }{2}{\upsilon }_{z,0}^2+\rho g{z}_0+\frac{\text{γ}}{r_0}-\frac{\text{γ}}{r\left(z\right)}=\frac{\rho }{2}{\upsilon }_z^2\left(z\right)+\rho gz\left(z\right)\end{array}$

$\begin{array}{ll}\frac{\rho }{2}{\upsilon }_{z,0}^2+\rho g{z}_0+\frac{\gamma }{r_0}\left(1-\frac{r_0}{r\left(z\right)}\right)\hfill & =\frac{\rho }{2}{\upsilon }_z^2\left(z\right)+\rho gz\left(z\right)\hfill \\ \frac{\rho }{2}{\upsilon }_{z,0}^2-\frac{\rho }{2}{\upsilon }_z^2\left(z\right)+\frac{\gamma }{r_0}\left(1-\frac{r_0}{r\left(z\right)}\right)\hfill & =\rho gz\left(z\right)-\rho g{z}_0\hfill \end{array}$

where we can set z = z − z₀ as the axis is chosen accordingly. We can further rewrite Eq. 23.7 to

$\begin{array}{ll}{\left(\frac{{\upsilon }_z^2\left(z\right)}{{\upsilon }_{z,0}}\right)}^2-1\hfill & =\frac{2}{\rho {\upsilon }_{z,0}^2}\left(gz+\frac{\gamma }{r_0}\left(1-\frac{r_0}{r\left(z\right)}\right)\right)\hfill \\ \frac{{\upsilon }_z^2\left(z\right)}{{\upsilon }_{z,0}}\hfill & ={\left(1+\frac{2}{\rho {\upsilon }_{z,0}^2}\left(gz+\frac{\gamma }{r_0}\left(1-\frac{r_0}{r\left(z\right)}\right)\right)\right)}^{\frac{1}{2}}\hfill \\ \frac{{\upsilon }_z^2\left(z\right)}{{\upsilon }_{z,0}}\hfill & ={\left(1+\frac{2g{r}_0}{{\upsilon }_{z,0}^2}\frac{z}{r_0}+\frac{2\gamma }{\rho {\upsilon }_{z,0}^2{r}_0}\left(1-\frac{r_0}{r\left(z\right)}\right)\right)}^{\frac{1}{2}}\hfill \\ \frac{{\upsilon }_z^2\left(z\right)}{{\upsilon }_{z,0}}\hfill & ={\left(1+\frac{2}{\text{Fr}}\frac{z}{r_0}+\frac{2}{\text{We}}\left(1-\frac{r_0}{r\left(z\right)}\right)\right)}^{\frac{1}{2}}\hfill \end{array}$

where we have used the Froude (see section 9.9.11) and the Weber (see section 9.9.12) numbers

$\begin{array}{ll}\text{Fr}\hfill & =\frac{{\upsilon }_{z,0}^2}{g{r}_0}\hfill \\ \text{We}\hfill & =\frac{\rho {\upsilon }_{z,0}^2{r}_0}{\gamma }\hfill \end{array}$

The Froude number allows assessing whether or not the effects of gravity can be neglected. This is the case for fluid jets with very high initial velocity. As we will see, these are hardly accelerated by gravitational forces. In these cases, the second term on the right-hand side of Eq. 23.8 can be set to 0. On the other hand, the Weber number allows assessing whether or not the effects of surface tension can be neglected. Again, this is the case for liquid jets with very high velocities for which the third term on the right-hand side of Eq. 23.8 can be set to 0.

Eq. 23.8 still contains the jet’s velocity as a function both of the independent variable z and the dependent variable r (z). Using the continuity equation (Eq. 10.7) we can establish a connection between the two:

$\begin{array}{ll}\dot{V}\hfill & =\pi r_0^2{\upsilon }_{z,0}={\displaystyle {\int }_0^{r\left(z\right)}{\upsilon }_z^2\left(z\right)2\pi rdr}\hfill \\ \hfill & ={\upsilon }_z^2\left(z\right)2\pi r^2\left(z\right)\hfill \\ \frac{{\upsilon }_z^2\left(z\right)}{{\upsilon }_{z,0}}\hfill & ={\left(\frac{r_0}{r\left(z\right)}\right)}^2\hfill \end{array}$

Inserting Eq. 23.9 into Eq. 23.8 results in

$\begin{array}{ll}{\left(\frac{r_0}{r\left(z\right)}\right)}^2\hfill & ={\left(1+\frac{2}{\text{Fr}}\frac{z}{r_0}+\frac{2}{\text{We}}\left(1-\frac{r_0}{r\left(z\right)}\right)\right)}^{\frac{1}{2}}\hfill \\ {\left(\frac{r_0}{r\left(z\right)}\right)}^2\hfill & 1+\frac{2}{\text{Fr}}\frac{z}{r_0}+\frac{2}{\text{We}}\left(1-\frac{r_0}{r\left(z\right)}\right)\hfill \end{array}$

Eq. 23.10 can be solved to find r (z). The result can be used to derive ${\upsilon }_z^2\left(z\right)$ using Eq. 23.9. However, Eq. 23.10 is not trivial to solve. The best method of obtaining a solution for it is by numerically solving for different values of z. Fig. 23.2 shows the result of a calculation on a jet of water that is pressed out of an opening with a diameter of 2 cm with a input overpressure p₀. This overpressure is varied from 10 bar to 1 mbar. Using Eq. 14.11 this pressure is converted to v_z,0 according to

$\begin{array}{ll}{\upsilon }_{z,0}\hfill & =\sqrt{\frac{2\left(p_0-p_{\infty }\right)}{\rho }}\hfill \\ \hfill & =\sqrt{\frac{2\Delta p}{\rho }}\hfill \end{array}$

with p_∞ = 1 bar and neglecting any differences in height across the exit of the nozzle. As can be seen from Fig. 23.2 the contraction of the fluid jet is strongly dependent on v_z,0 and therefore on ∆p. ∆p = 10 bar is the maximum pressure that is provided by the water tap connected to the drinking water supply. In this case, Fr ≈ 20 000 and We ≈ 250 000, and therefore the inertia terms strongly dominate over both gravity and surface tension forces. In consequence, there is virtually no contraction of the water jet, as can be seen from Fig. 23.2. The right-hand side of Eq. 23.9 reduces to only the first term as both Fr and We are high.

On reducing ∆p the influence of the gravitation force becomes visible. For ∆p = 100 mbar, Fr ≈ 2000 and We ≈ 27 000. The effects due to surface tensions are still one order of magnitude smaller than the effects of gravitational forces. The influence of the latter can be seen in the linear decrease of the diameter with increasing z (see Eq. 23.10). The jet will contract to about 95 % of its original diameter over a total distance of 1000 mm.

For ∆p = 100 mbar, the increasing radius contraction becomes obvious. Over a distance of 1000 mm, the radius reduces to roughly 50 % of its original radius. In this case Fr ≈ 200 and We ≈ 2700. Here, the effects due to surface tension (We) are acting in addition to the effects of gravitation, which result in a nonlinear contraction profile along z.

For ∆p = 10 mbar the jet radius after 1000 mm is only about 10 % of the original radius. The corresponding dimensionless numbers are Fr ≈ 20 and We ≈ 270. Here the effect of both gravitation and surface tension is strongly influential on the fluid mechanics of the jet.

Finally, for ∆p = 1 mbar the contraction will reduce the jet diameter to only a few percent of its original size. In this case, Fr ≈ 2 and We ≈ 27. How one can clearly see that both gravitational and surface tension forces are in the same orders of magnitude as the inertia forces. In this regime, the jet will not be stabilized by its inertia and will eventually decompose as a consequence of the gravitational acceleration and the surface tension forces.

23.4 Instability

Derivation. Due to the continuous decrease of the radius, fluid jets will eventually become instable. This effect is not only monitored on jets accelerated by gravity, but also in closed microfluidic systems where a flow is squeezed in between two immiscible sheathing streams as used, e.g., in droplet microfluidics for the generation of droplets. In both cases, the fluid jet will eventually become instable and breakup into individual droplets. As we have seen, surface tension forces play a significant role in this effect.

Fig. 23.3a details the process of the jet instability. After a given distance the jet will show small ripplets that originate from oscillations, e.g., in the driving pressure. If these oscillations are at the correct frequency they will augment in amplitude and result in the development of visible perturbations that will eventually destabilize the jet. The result is the formation of droplets that are, as discussed in section 23.2, thermodynamically more stable. As noted, these effects do not only occur on the falling fluid jet. As we will see, this destabilization is mainly due to the increase in Young-Laplace pressure gradients induced by the formation of curved surfaces on the jet. These curvatures are the result of the perturbation acting on the fluid jet. These perturbations are usually only very minor in size, i.e., only very tiny changes in the pressure, purity, or composition of the jet’s fluid will be sufficient to induce small perturbations. Usually these perturbations will decay quickly. However, if they occur at a given frequency, they will amplify similar to the formation of a resonance in solids when excited at the resonance frequency.

Analytically, we assume the falling jet as discussed in section 23.3. Its radius r (z) is dependent on the independent variable z, as is its radial velocity v_r (z) and its axial velocity v_z (z), which we already discussed in section 23.3. The last variable of note is the pressure inside of the jet p (z), which is also dependent on z. We will now assume a cyclic perturbation on the radius of the jet diameter, which is expressed as

$\begin{array}{ll}r\left(z,t\right)\hfill & =r_0\in {\text{e}}^{\omega t}{\text{e}}^{\text{i}kz}\hfill \\ \hfill & =r_0\in {\text{e}}^{\omega t}{{}^{\text{+}}}^{\text{i}kz}\hfill \end{array}$

The perturbation has a radial component of amplitude r₀ that travels (in time) at the angular frequency ω and along the z-axis at a speed given by the wave number k. As these two velocity components are orthogonal, i.e., in radial and axial directions, we formulate the radial disturbance as a complex value. As we have seen, the radius along the z-axis will change (see Eq. 23.10). However, as the wavelength of the disturbance is very small, we can approximate the radius of the unperturbed jet r₀ around a given position z to be approximately constant.

The perturbation in radius will also lead to perturbations in the other three variables that should occur at the same frequency, and therefore

$r\left(z,t\right)=r_0+\in {\text{e}}^{\omega t}{{}^{\text{+}}}^{ikz}$

${\upsilon }_r\left(r,z,t\right)={\upsilon }_{r,0}\left(r\right){\text{e}}^{\omega t}{{}^{\text{+}}}^{ikz}$

${\upsilon }_z\left(r,z,t\right)={\upsilon }_{z,0}\left(r\right){\text{e}}^{\omega t}{{}^{\text{+}}}^{ikz}$

$p\left(r,z,t\right)=p_0+p_0\left(r\right){\text{e}}^{\omega t}{{}^{\text{+}}}^{ikz}$

In all cases we assume only minor changes around a given z value in which case we can assume the unperturbed parameter to be approximately constant. We now apply the Navier-Stokes equation in cylindrical coordinates assuming all azimuthal changes to be zero since we only introduced radial perturbations. In this case, the Navier-Stokes equation (see Eq. 13.2) reads

$\frac{\partial \overrightarrow{\upsilon }}{\partial t}+\upsilon \frac{\partial \overrightarrow{\upsilon }}{\partial \overrightarrow{s}}=-\frac{1}{\rho }\frac{\partial p}{\partial \overrightarrow{s}}$

in radial direction

$\frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial t}+{\upsilon }_r\left(r,z,t\right)\frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial r}+{\upsilon }_z\left(r,z,t\right)\frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial r}=-\frac{1}{\rho }\frac{\partial p}{\partial r}$

in axial (z) direction

$\frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial t}+{\upsilon }_r\left(r,z,t\right)\frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial z}+{\upsilon }_z\left(r,z,t\right)\frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}$

for which we can find the partial derivatives from Eq. 23.11 to Eq. 23.14 as

$\begin{array}{l}\frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial t}=\omega {\upsilon }_{r,0}\left(r\right)e^{\omega t+ikz}\\ \frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial t}=\omega {\upsilon }_{z,0}\left(r\right)e^{\omega t+ikz}\\ \frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial r}=\frac{\partial {\upsilon }_{r,0}\left(r\right)}{\partial r}{e}^{\omega t+ikz}=\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}{e}^{\omega t+ikz}\\ \frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial r}=\frac{\partial {\upsilon }_{z,0}\left(r\right)}{\partial r}{e}^{\omega t+ikz}=\frac{d{\upsilon }_{z,0}\left(r\right)}{dr}{e}^{\omega t+ikz}\\ \frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial z}=ik{\upsilon }_{r,0}\left(r\right)e^{\omega t+ikz}\\ \frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial r}=ik{\upsilon }_{z,0}\left(r\right)e^{\omega t+ikz}\\ \frac{\partial p\left(r,z,t\right)}{\partial r}=\frac{\partial p_0\left(r\right)}{\partial r}{e}^{\omega t+ikz}=\frac{\partial p_0\left(r\right)}{dr}{e}^{\omega t+ikz}\\ \frac{\partial p\left(r,z,t\right)}{\partial z}=ik{p}_0\left(r\right)e^{\omega t+ikz}\end{array}$

which we can substitute into Eq. 23.15 and Eq. 23.16, resulting in the following:

in radial direction

${\text{e}}^{\omega t+\text{i}kz}\left({\upsilon }_{r,0}\left(r\right)\omega +{\upsilon }_{r,0}\left(r\right)\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}{\text{e}}^{\omega t+\text{i}kz}+{\upsilon }_{z,0}\left(r\right)\frac{d{\upsilon }_{z,0}\left(r\right)}{dr}{\text{e}}^{\omega t+\text{i}kz}\right)=-\frac{1}{\rho }\frac{d{p}_0\left(r\right)}{dr}{\text{e}}^{\omega t+\text{i}kz}$

in axial (z) direction

${\text{e}}^{\omega t+\text{i}kz}\left({\upsilon }_{z,0}\left(r\right)\omega +\text{i}k{\upsilon }_{r,0}{\left(r\right)}^2{\text{e}}^{\omega t+\text{i}kz}+\text{i}k{\upsilon }_{z,0}{\left(r\right)}^2{\text{e}}^{\omega t+\text{i}kz}\right)=-\frac{1}{\rho }ik{p}_0\left(r\right){\text{e}}^{\omega t+\text{i}kz}$

If we assume the perturbation to be small, we can simplify Eq. 23.17 and Eq. 23.18 by neglecting all terms that contain (e^{ωt+i k z})² which results in the following:

in radial direction

${\upsilon }_{r,0}\left(r\right)\omega =-\frac{1}{\rho }\frac{d_{p0}\left(r\right)}{dr}$

in axial (z) direction

${\upsilon }_{z,0}\left(r\right)\omega =-\frac{1}{\rho }\text{i}k{p}_0\left(r\right)$

applying $\frac{d}{dr}$ on both sides

$\begin{array}{ll}\frac{d}{dr}\left({\upsilon }_{z,0}\left(r\right)\omega \right)\hfill & =\frac{d}{dr}\left(-\frac{1}{\rho }\text{i}k{p}_0\left(r\right)\right)\hfill \\ \frac{d{\upsilon }_{z,0}\left(r\right)}{dr}\omega \hfill & =-\frac{1}{\rho }\text{i}k\frac{d{p}_0\left(r\right)}{dr}\hfill \\ \frac{d{p}_0\left(r\right)}{dr}\hfill & =-\frac{\rho \omega }{\text{i}k}\frac{d{\upsilon }_{z,0}\left(r\right)}{dr}\hfill \end{array}$

We can now insert Eq. 23.20 into Eq. 23.19 and obtain

$\begin{array}{ll}{\upsilon }_{r,0}\left(r\right)\hfill & =\frac{1}{\text{i}k}\frac{d{u}_{z,0}\left(r\right)}{dr}\hfill \\ \frac{d{\upsilon }_{z,0}\left(r\right)}{dr}\hfill & =\text{i}k{\upsilon }_{r,0}\left(r\right)\hfill \end{array}$

We have now found an analytical correlation between v_r,0 (r) and v_z,0 (r). Introducing the continuity equation in cylindrical form (see Eq. 13.1) allows reducing the partial differential equation (Eq. 23.19) to an ordinary differential equation. Again, we neglect any changes in the azimuthal direction as the perturbation is only in the radial direction. The continuity equation then simplifies to

$\frac{\partial {\upsilon }_r\left(r,z,t\right)}{\partial r}+\frac{{\upsilon }_r\left(r,z,t\right)}{r}+\frac{\partial {\upsilon }_z\left(r,z,t\right)}{\partial z}=0$

into which we can insert the partial differential that we formulated earlier to result in

$\begin{array}{ll}\frac{\partial {\upsilon }_{r,0}\left(r\right)}{\partial r}{\text{e}}^{\omega t+\text{i}kz}+\frac{{\upsilon }_{r,0}\left(r\right){\text{e}}^{\omega t+\text{i}kz}}{r}+\text{i}k{\upsilon }_{z,0}\left(r\right){\text{e}}^{\omega t+\text{i}kz}\hfill & =0\hfill \\ \frac{\partial {\upsilon }_{r,0}\left(r\right)}{\partial r}+\frac{{\upsilon }_{r,0}\left(r\right)}{r}+\text{i}k{\upsilon }_{z,0}\left(r\right)\hfill & =0\hfill \\ \frac{\partial }{\partial r}\left(\frac{\partial {\upsilon }_{r,0}\left(r\right)}{\partial r}+\frac{{\upsilon }_{r,0}\left(r\right)}{r}+\text{i}k{\upsilon }_{z,0}\left(r\right)\right)\hfill & =0\hfill \\ \frac{d^2{\upsilon }_{r,0}\left(r\right)}{d{r}^2}-\frac{{\upsilon }_{r,0}\left(r\right)}{r^2}+\frac{1}{r}\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}+\text{i}k\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}\hfill & =0\hfill \end{array}$

Inserting Eq. 23.21 into Eq. 23.22 yields the second-order homogeneous ODE:

$\begin{array}{l}\frac{d^2{\upsilon }_{r,0}\left(r\right)}{d{r}^2}-\frac{{\upsilon }_{r,0}\left(r\right)}{r^2}+\frac{1}{r}\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}-k^2{\upsilon }_{r,0}\left(r\right)=0\hfill \\ r^2\frac{d^2{\upsilon }_{r,0}\left(r\right)}{d{r}^2}+r\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}-{\upsilon }_{r,0}\left(r\right)-r^2{k}^2{\upsilon }_{r,0}\left(r\right)=0\hfill \\ r^2\frac{d^2{\upsilon }_{r,0}\left(r\right)}{d{r}^2}+r\frac{d{\upsilon }_{r,0}\left(r\right)}{dr}-\left(1+{\left(kr\right)}^2\right){\upsilon }_{r,0}\left(r\right)=0\hfill \end{array}$

Eq. 23.23 is a modified Bessel function (see section 3.2.4.2) of first order (ν = 1). It is solved by a function of type Eq. 3.69 resulting in

${\upsilon }_{r,0}\left(r\right)=c_1{I}_1\left(kr\right)$

Using this solution for v_r,0 (r) (for which we still need to find the integration constant c₁) allows deriving v_z,0 (r) from Eq. 23.21, which we can use to solve Eq. 23.20 for p₀ (r). As Eq. 23.23 is a second-order differential equation we will need a total of two boundary conditions.

Boundary Condition 1: Kinematic Boundary Condition. The first boundary condition is given by the radial speed of the fluid on the surface of the jet. This velocity is given only by the temporal change of the diameter of the jet. As there is no net flow across the boundary of the jet (otherwise the conservation of mass would not be fulfilled) the only velocity component is given by the velocity of the boundary itself. This boundary condition is therefore referred to as the kinematic boundary condition. Using Eq. 23.11 we can derive

$\begin{array}{ll}\frac{\partial }{\partial t}\left(r_0\left(z,t\right)\right)\hfill & =\frac{\partial }{\partial t}\left(r_0+\in {\text{e}}^{\omega t+\text{i}kz}\right)\hfill \\ \hfill & =\in \omega {\text{e}}^{\omega t+\text{i}kz}={\upsilon }_r\left(r_0\right){\text{e}}^{\omega t+\text{i}kz}\hfill \end{array}$

Using Eq. 23.24 and Eq. 23.25 we can now derive c₁ as

$\begin{array}{ll}{c}_1{I}_1\left(k{r}_0\right)\hfill & =\in \omega \hfill \\ c_1\hfill & =\frac{\in \omega }{I_1\left(k{r}_0\right)}\hfill \end{array}$

Substituting Eq. 23.20 into Eq. 23.21 we find

$\frac{d{p}_0\left(r\right)}{dr}=-\rho \omega {\upsilon }_{r,0}\left(r\right)$

Combining Eq. 23.25 with Eq. 23.27 we find

$\frac{d{p}_0\left(r\right)}{dr}=-\rho \omega {\upsilon }_{r,0}\left(r\right)$

into which we can insert Eq. 23.24

$\frac{d{p}_0\left(r\right)}{dr}=-\rho \omega c_1{I}_1\left(kr\right)$

We can now integrate Eq. 23.28 to result in

$\begin{array}{ll}{p}_0(r)\hfill & =-\frac{\rho \omega }{k}{c}_1{I}_0(kr)\hfill \\ \hfill & =-\frac{\rho \in {\omega }^2}{k}\frac{I_0(kr)}{I_1(k{r}_0)}\hfill \end{array}$

where we have used the recurrence relationship of the modified Bessel function (see Eq. 3.74) as well as the value for c₁ (see Eq. 23.26). We can now express the pressure relation using Eq. 23.14 as

$p\left(r,z,t\right)=p_0-\frac{\rho \in {\omega }^2}{k}\frac{I_0\left(kr\right)}{I_1\left(k{r}_0\right)}{\text{e}}^{\omega t+\text{i}kz}$

Boundary Condition 2: Pressure Gradient Induced by Surface Curvature. We can derive a second solution for p (r, z, t) using the Young-Laplace pressure gradient induced by the doubly curved surface. For this, we need the two radii, r₁ and r₂. Radius r₁ can be derived from as Eq. 23.11

$\frac{1}{r_1}=\frac{1}{r_0+\in {\text{e}}^{\omega t+\text{i}kz}}\approx \frac{1}{r_0}-\frac{\in {\text{e}}^{\omega t+\text{i}kz}}{r_0^2}$

where we applied a little calculus trick in order to simplify the denominator. For details on this simplification used see Eq. 3.95. We can use this trick as ω e^{ω t+i k z} r₀. Radius r₂ can be found by using the curvature κ (see Eq. 3.103)

$\begin{array}{ll}\frac{1}{r_2}=-{\kappa }_2\hfill & =-\frac{d^{2}r\left(z,t\right)}{d{z}^2}\hfill \\ \hfill & =k^2\in {\text{e}}^{\omega t+\text{i}kz}\hfill \end{array}$

Please note that r₁> 0 since it is a convex surface, whereas r₂< 0 since it is a concave radius. There is no reason for choosing whether the concave or the convex radii are positive; it is simply important to ensure that one of them is negative.

Using Eq. 20.11 we can now find the Young-Laplace pressure drop between the perturbed and the unperturbed fluid cylinder as

$p\left(r,z,t\right)=\gamma \left(\left(\frac{1}{r_{1,\text{perturbed}}}+\frac{1}{r_{1,\text{perturbed}}}\right)-\left(\frac{1}{r_{2,\text{unperturbed}}}+\frac{1}{r_{2,\text{unperturbed}}}\right)\right)+p_0$

with

$\begin{array}{l}{r}_{1,\text{unperturbed}}=r_0\hfill \\ r_{2,\text{unperturbed}}=\infty \hfill \end{array}$

which leads to

$\begin{array}{ll}p\left(r,z,t\right)\hfill & =\gamma \left(\frac{1}{r_{1,\text{perturbed}}}+\frac{1}{r_{2,\text{perturbed}}}-\frac{1}{r_0}\right)+p_0\hfill \\ \hfill & =\gamma \left(\frac{1}{r_0}-\frac{\in {\text{e}}^{\omega t+\text{i}kz}}{r_0^2}+k^2\in {\text{e}}^{\omega t+\text{i}kz}-\frac{1}{r_0}\right)+p_0\hfill \\ \hfill & =\frac{\gamma \in }{r_0^2}\left(1-k^2{r}_0^2\right)\in {\text{e}}^{\omega t+\text{i}kz}+p_0\hfill \end{array}$

With Eq. 23.30 we have found a second relationship for the pressure inside the perturbed fluid jet.

Dispersion Relation. In order to balance the forces inside of the fluid, the pressure induced by the perturbation must be balanced by the surface tension forces. Therefore, Eq. 23.29 must be equal to Eq. 23.30. This results in the dispersion relation. Dispersion relations are equations that describe the dependency of a physical system as a function of the variables. Here, the dispersion relation is given by

$\begin{array}{ll}{p}_0-\frac{\rho \in {\omega }^2}{k}\frac{I_0\left(kr\right)}{I_1\left(k{r}_0\right)}{\text{e}}^{\omega t+\text{}\text{i}kz}\hfill & =p_0-\frac{\gamma \in }{r_0^2}\left(1-k^2{r}_0^2\right){\text{e}}^{\omega t+\text{}\text{i}kz}+p_0\hfill \\ \frac{\rho \in {\omega }^2}{k}\frac{I_0\left(kr\right)}{I_1\left(k{r}_0\right)}\hfill & =\frac{\gamma \in }{r_0^2}\left(1-k^2{r}_0^2\right)\hfill \\ \frac{\rho r_0^3}{\gamma }{\omega }^2\hfill & =k{r}_0\left(1-k^2{r}_0^2\right)\frac{I_1\left(k{r}_0\right)}{I_0\left(kr\right)}\hfill \end{array}$

The dispersion relation, Eq. 23.31, allows us to study at which conditions the system will become instable. If the right-hand side of the equation becomes large, ω has to become large as the term $\frac{\rho r_0^3}{\text{γ}}$ is constant. This is $k^2{r}_0^2=1$, which is only possible if

$\begin{array}{ll}k{r}_0\hfill & <1\hfill \\ \frac{2\pi }{\lambda }{r}_0\hfill & <1\hfill \\ \lambda \hfill & =2\pi r_0\hfill \end{array}$

Therefore, the system can become instable if the wavelength of the perturbation is bigger than the diameter of the jet. As we have seen from Eq. 23.10, the radius of a falling jet will decrease continuously, in which case the instability condition, Eq. 23.32, will eventually be fulfilled.

Fig. 23.4 shows a plot of the dispersion relation given by Eq. 23.32. The function reaches its maximum at k r = 0.697, for which the value reached is $\frac{\rho r_0^3}{\text{γ}}{\omega }^2\approx 0.1178$. Therefore the fastest growing perturbation is the perturbation with a wavelength of $\lambda =\frac{2\pi }{0.697}{r}_0=9.01r_0$.

23.5 Standing Waves on a Fluid Jet

The dispersion relation (Eq. 23.31) allows explaining another commonly encountered phenomenon. When a fluid jet drops into a large reservoir or onto a solid stationary surface, e.g., a finger, one can observe a small zone of seemingly standing waves excited some millimeters away from the stationary surface. These waves are actually perturbations that travel at the exact negative speed of the fluid jet. Thus they appear, to the stationary observer, to be stationary in space (see Fig. 23.5). In this case

${\upsilon }_{\text{jet}}=-\frac{\omega }{k}$

If v_jet is known, inserting Eq. 23.33 into Eq. 23.32 will allow solving for the wavelengths that will appear to be stationary. If v_jet is not known, it can be determined by measuring the wavelengths of the stationary perturbations and inserting this value into Eq. 23.32. This is a commonly used principle for measurements of the jet’s speed.

The position of the stationary perturbations is a function of the fluid’s surface tension. It can be altered by dropping surfactants into the fluid reservoir or by wetting the finger with soap. Again, this phenomenon can be derived directly from Eq. 23.32.

23.6 Characteristic Breakup Time

The characteristic breakup time tbreakup,char. of the jet can be found as the inverse of the maximum angular frequency ω_max for k r = 0.697:

$\begin{array}{ll}\frac{\rho r_0^3}{\gamma }{\omega }^2\hfill & 0.1178\hfill \\ t_{\text{breakup,}\text{char}.}\hfill & =\frac{1}{\omega }=\sqrt{0.1178\frac{\gamma }{\rho r_0^3}}\hfill \\ \hfill & =2.91\sqrt{\frac{\rho r_0^3}{\gamma }}\hfill \end{array}$

Fig. 23.6 displays the characteristic breakup time of fluid jets as a function of the jet’s radius. For a jet of diameter 1 cm, the breakup time is approximately 0.121 s. Obviously, the critical length can also be calculated if the average velocity of the jet is known.

23.7 Applicability of the Plateau-Rayleigh Instability

There are a number of practical applications for the theory of the Plateau-Rayleigh instability phenomenon. Liquid dispensers, inkjet printers, and similar spotting systems are designed to induced perturbations (either by means of mechanical oscillations or thermally) in order to ensure correct decomposition of a fluid jet into droplets. As an everyday phenomenon, consider the fact that rain does not fall as thin cylinders of water but rather in droplet form. Obviously, droplet microfluidics makes extensive use of this concept for the creation of small droplets. Another important application scope for the instability is the breakup of thin layers of viscous fluids, e.g., glues or paints on planar surfaces and wires. The layers tend to break up during retraction of the surface from the liquid, which results in inhomogeneous coatings.

One important thing to keep in mind is that, during the fluid mechanical derivation of the instability, several restrictions have been made. The most important one is that the viscous effects in the fluid have largely been neglected compared to the inertia and surface tension forces. During the derivation of the falling jet (see section 23.3) we have already introduced the Froude and the Weber numbers that allowed assessing the influence of inertia forces compared to gravitational (Froude) and surface tension forces (Weber). During this derivation we neglected viscous forces. The dimensionless number that indicates whether or not this is appropriate is the Ohnesorge number (see section 9.9.13), which puts viscous forces in relation to both inertia and surface tension forces.

Using water, the range of the Ohnesorge number is very small if the characteristic length is above roughly 1 µm (see Fig. 23.7). However, for different liquids, e.g., inkjet printer inks, that are mixtures glycols, e.g., ethylene glycol and aqueous inks, the viscosity values are higher (usually somewhere above 7 mPa s) and the surface tension values are lower (usually in the range of about 40 mN m⁻¹). Therefore, the Ohnesorge numbers may become larger and thus the experimental values may deviate from the theoretical values (see Fig. 23.7).

23.8 Summary

In this section, we have studied the Plateau-Rayleigh instability, which is the mechanism by which initially compact liquid jets will eventually rupture and decompose into droplets. This phenomenon is commonly observed in everyday life, e.g., when observing water running from a faucet. As we have seen, the characteristic breakup time depends on several factors, including the diameter of the jet and the velocity at which the liquid propagates.