We examine the flow in and around a falling fluid droplet in a vertically oscillating flow. We assume axisymmetric Stokes flow, and for small deformations to the droplet, the governing equations can be linearized leading to an infinite system of linear ordinary differential equations. In this study, we have analytically solved the problem in the small-capillary limit. We note that the solution locally breaks down at the poles of the droplet. The drag and center of the mass were also obtained. In the case when only odd modes are present, the droplet shows three-dimensional axisymmetric heart-shaped solutions oscillating vertically in time. When only even modes are present, the droplet exhibits axisymmetric stretching and squeezing.

The classical problem of a sphere moving through a fluid has long attracted great interest for its numerous applications but also fundamental significance for fluid dynamics. For example, Chatterjee et al.1 examined whether the flow from an expanding or collapsing microbubble near a cell could be used as a drug delivery technique, Ward et al.2 considered whether such flows could be used to destroy cancerous cells, and Krehbiel et al.3 investigated whether they could be used to rupture algal cells. Other applications include sedimentation, lubrication processes, emulsions, and suspensions, for example, microorganisms, paint, sun protection cream, etc. As the problem has a long history, we shall only highlight some of the key works.

Almost two centuries ago, Stokes4 examined steady flow past a solid sphere of radius Rc and moving at uniform speed U in the absence of inertia and he obtained the stream function of the flow field $ψ = ( U R c / 4 r ) ( 3 r 2 − R c 2 ) sin 2 ( θ )$, where r is the radial coordinate measured from the center of the sphere and θ is the angle measured from the axis in the direction of the flow. Furthermore, he found that the drag on the sphere in the vertical direction was given by $D f = 6 π μ R c U$, where μ is the dynamic viscosity of the fluid. Stokes4 also solved the problem of a rigid sphere oscillating within a fluid in a spherical container in terms of the streamfunction using separation of variables. The so-called “Stokes solution” is one of the fundamental results in low-Reynolds-number hydrodynamics.

Oseen5 obtained the first-order correction to flow past a solid sphere for low Reynolds numbers, with the dimensionless streamfunction given by
$ψ = 2 r 3 + 1 4 r sin 2 ( θ ) − 3 ( 1 + cos ( θ ) ) 2 Re ( 1 − e − r 2 Re ( 1 − cos ( θ ) ) ),$
and, further, that the drag coefficient on the sphere in the vertical direction was $C d = 6 π Re ( 1 + 3 Re 8 )$. Proudman and Pearson6 obtained a higher-order approximation to flow past a solid sphere for low Reynolds numbers, Re, which includes corrective terms of $O ( Re 2 ln ( Re ) )$, and they found the drag coefficient was given by $C d = 6 π Re ( 1 + 3 Re 8 + 9 40 Re 2 ln ( Re ) + O ( Re 2 ) )$. Payne and Pell7 explored Stokes flow for a class of axially symmetric solid bodies and obtained the drag on a variety of bodies including a lens-shaped body, hemisphere, spherical cap, a pair of separated sphere, a spheroid, and a lens. Cox8 obtained the drag in the low Reynolds number limit up to $O ( Re 2 ln ( Re ) )$ for steady flow around arbitrary-shaped solid bodies falling at a constant speed, such bodies included a moving spheroid, a moving dumb-bell-shaped body, a moving rotating sphere, and a dumb-bell-shaped body in pure rotation. Ockendon9 considered unsteady flow past a solid sphere with a time-dependent velocity at small-but-finite Reynolds numbers and showed that the drag predicted by the Stokes flow differs from that obtained from the unsteady Navier–Stokes solution. Chester et al.10 considered flow past a solid sphere for low Reynolds numbers, Re, which includes corrective terms of $O ( Re 3 ln ( Re ) )$, and they found the drag coefficient was given by $C d = 6 π Re ( 1 + 3 Re 8 + 9 40 Re 2 [ ln ( Re ) + c ] + 27 80 Re 3 ln ( Re ) + O ( Re 2 ) )$, where $c = γ + 5 3 ln ( 2 ) − 323 360$ and γ is Euler's constant. On the other hand, Pruppacher et al.11 numerically examined flow past a solid sphere for moderate Reynolds numbers, which agreed well with experimentally obtained values of the drag. The numerical results for the drag agreed well with the analytical results obtained in the low Reynolds number limit, for small Reynolds numbers, but these analytical results diverged from the numerical solution for moderate Reynolds numbers.

Landau and Lifshitz12 gave the solution for an oscillating spherical drop in an infinite medium and found the smallest possible frequency of oscillations of the drop was $8 α / ( ρ R c 3 )$, where α is the surface tension coefficient, ρ is the density of the fluid, and Rc is the radius of the droplet. They said that “the oscillations cause the surface of the drop to deviate from the spherical form.” Mei and Adrian13 examined unsteady low Reynolds number with very low-frequency oscillatory flow past a stationary solid sphere and found that the acceleration-dependent force was linearly proportional to the frequency. They found that the classical Stokes solution was not valid for small frequencies for small Reynolds numbers. Chang et al.14 focused on axisymmetric viscous laminar flow around solid spheroids for moderate Reynolds numbers. They found that, for small times, the asymptotic analysis and numerical solutions obtained using finite differencing agreed well. Taseli and Demiralp15 examined axisymmetric Stokes flow past an arbitrary axisymmetrical solid body by writing the solution as an infinite series involving Gegenbauer polynomials. Otto16 explored the stability of the flow around a solid sphere oscillating at a high frequency. The problem was reduced to an infinite system of ordinary differential equations. Using linear stability analysis, they found that the flow could become unstable to Taylor–Görtler vortices.

The flow of an fluid sphere through a another fluid in the absence of inertia was analyzed by Rybczynski17 and Hadamard18 who independently found that the dimensionless streamfunction inside $ψ ̂$ and outside ψ of the sphere were given by
$ψ ̂ = 3 − 2 B 4 ( r 4 − r 2 ) sin 2 ( θ ) , ψ = 1 2 ( r 2 − B r + B − 1 r ) sin 2 ( θ ) ,$
where $B = ( 2 μ + 3 μ ̂ ) / ( 2 μ + 2 μ ̂ )$, where $μ ̂$ and μ denote the dynamic viscosity's of the fluids inside and outside the sphere, respectively. They found that the drag on the sphere was given by $D f = 4 3 π ( ρ ̂ − ρ ) g R c 3$ where Rc is the radius of the sphere, g is the magnitude of the gravitational acceleration, and $ρ ̂$ and ρ denote the density of the fluids inside and outside the sphere, respectively. Furthermore, the speed of the sphere was given by $U = g R c 2 ( ρ ̂ − ρ ) / ( 3 B μ )$. Taylor and Acrivos19 theoretically investigated the axisymmetric motion of a slightly deformable fluid drop falling through a fluid in the small-but-finite Reynolds number limit. They found that for small Weber numbers, the drop will deform into an oblate spheroid while further increase in the Weber number deforms the droplet into a spherical cap shape. Lin and Gautesen20 studied the small-but-finite Reynolds number flow of axisymmetric steady fluid surrounding a deformable sphere with variable radius. They obtained the drag up to $O ( Re 2 ln ( Re ) )$. To illustrate there result by considering two cases: a pulsating sphere and a constantly expanding sphere. Oliver and Chung21 numerically considered flow inside and outside a fluid sphere at low Reynolds number for a variety of density ratios. They found that the drag increases when the viscosity ratio is increased, but decreased when the Reynolds number was increased. They found that the density ratio had little effect on the drag. Pozrikidis22 examined a viscous drop subject to axisymmetric perturbations. They found that a moving spherical drop was unstable and developed into a nearly steady ring under perturbations. Furthermore, surface tension was not capable of suppressing the instability. Machu et al.23 numerically and experimentally examined the small-but-finite Reynolds number flow around a deforming droplet. They found that everything they observed experimentally could be observed using Stokes flow without the need to include surface tension or inertial effects. Srivastava et al.24 numerically investigated the steady flow around an oblate axisymmetric body for various eccentricities. They found that increasing the eccentricity of the deformed sphere reduced the drag with a flat circular disk having the smallest drag. Krehbiel and Freund25 considered axisymmetric steady inviscid flow surrounding a Newtonian liquid sphere. They were able to obtain analytical solutions for the inner and outer streamfunction as relatively simple finite expressions. Recently, Sahu and Khair26 numerically investigated a neutrally buoyant viscous droplet and found that the droplet could break up if the capillary number was greater than a critical value that depended on the Deborah number. Furthermore, Godé et al.27 numerically examined the Basset–Boussinesq history force on a droplet in a uniform oscillating flow. By adjusting the frequency of the oscillation, they were able to determine the range of physical parameters that make the contribution of the history force significant.

There are several studies involving non-Newtonian fluids. Leslie and Tanner28 examined low-Reynolds-number flow of an axisymmetric steady non-Newtonian fluid surrounding a solid sphere. They found the drag on a solid sphere by the non-Newtonian fluid to be smaller than the drag on a solid sphere by a Newtonian fluid. Caswell and Schwarz29 looked at low-Reynolds-number flow of a non-deformable Newtonian spherical droplet surrounded by an incompressible Rivlin–Ericksen fluid. Sadly, their analytical expression for the drag on the sphere involved two unknown parameters, which could not be obtained from the experimental data available to them, so they were unable to compare their work with previous studies. Beris et al.30 considered a solid sphere falling through a Bingham plastic material. They numerically solved the flow field using the finite element method. They found that the drag on the sphere was greater in a Bingham plastic material, compared to than the drag on a solid sphere by a Newtonian fluid. By obtaining the drag on a sphere, one may be able to determine various physical properties of the fluid. Ramkissoon31 analytically examined steady axisymmetric Stokes flow past a non-deformable Reiner–Rivlin fluid spheroid. They obtained an analytical expression for the drag when the spheroid is only a slightly deformed sphere. They found the drag on the Reiner–Rivlin fluid spheroid is less than the drag on a Newtonian spheroid. Sostarecz and Belmonte32 experimentally examined an Order Three (see Bird et al.33) non-Newtonian fluid droplet falling through a Newtonian fluid. The droplet was found to exhibit a stable dimple at its edge, with the dimple moving toward the center of the droplet as the droplet volume increases, eventually leading to a torus-shaped droplet for sufficiently large droplet volumes. Mukherjee and Sarkar34 numerically investigated the motion of an Oldroyd-B fluid droplet falling in a Newtonian fluid using finite differences. They found the flow to be unstable when there was a decrease in surface tension. Jaiswal and Gupta35 analytically examined axisymmetric steady Stokes flow surrounding a Reiner–Rivlin liquid spheroid, which is very close to a sphere in shape. They obtained the flow field and drag on the spheroid. They found that the drag on a solid spheroid is greater than the drag on a Reiner–Rivlin liquid spheroid. Furthermore, the drag on a Reiner–Rivlin liquid spheroid is greater than the drag on a liquid sphere.

Vamerzani et al.36 analytically examined a deformable fluid droplet falling through a fluid using Stokes flow. They found good agreement between analytical and experimental results in estimating the terminal velocity and drop shape when both the Deborah and capillary numbers were small. Interestingly, it was observed that as the volume of the drop increases, the drop loses its spherical shape and falls faster. We note that some of the cross sections of their droplets resemble heart shapes. Norouzi and Davoodi37 investigated slightly deformable spherical droplets in Stokes flow when both the Deborah and capillary numbers were small. Again some of the droplets resemble heart shapes. The results were compared with experiments involving a fluid droplet falling through a fluid when both fluids were Oldroyd-B fluids. Jaiswala38 explored the axisymmetric steady motion of a Reiner–Rivlin fluid surrounding a Newtonian liquid spheroid, which is very close to a sphere in shape. For fluids with a smaller viscosity ratio, the droplet's speed will initially increase and then decrease as a function of the Weissenberg number.

In the present study, we examine axisymmetric Stokes flow in and around a falling fluid droplet under external forcing. In Sec. II, we present the problem and non-dimensionalize the governing equations and boundary conditions. In Sec. III, the equations are expressed in axisymmetric spherical polar coordinate while also introducing appropriate streamfunctions, and are linearized assuming the droplet is only slightly deformable. Section IV gives the well-known non-deformable droplet solutions. An infinite system of equations that the first-order (in terms of the droplet deformation parameter) solutions need to satisfy is derived in Sec. V. This system is rescaled in Sec. VI. In Sec. VII, expressions for the drag on the droplet in the vertical direction, the volume, and center of mass of the droplet are obtained. In Sec. VIII, we obtain the first-order steady-state solution in the small-capillary limit. In Sec. IX, we obtain the first-order unsteady solution in the small-capillary limit. Finally, a summary of our findings and conclusions are offered in Sec. X.

Suppose we have a droplet of fluid 2 falling through an infinite region of fluid 1 as illustrated in Fig. 1. We suppose that both are immiscible, incompressible Newtonian fluids of constant density and constant kinematic viscosity. The bulk equations are the Navier–Stokes and continuity equations
$∂ U ∂ T + ( U · ∇ ) U = − 1 ρ ∇ P + ν ∇ 2 U − g e z , ∇ · U = 0,$
(1)
$∂ U ̂ ∂ T + ( U ̂ · ∇ ) U ̂ = − 1 ρ ̂ ∇ P ̂ + ν ̂ ∇ 2 U ̂ − g e z , ∇ · U ̂ = 0,$
(2)
where hats denote fluid 2 (the droplet), U is the fluid velocity, T is time, ρ is the density, P is the pressure, ν is the kinematic viscosity, g is the magnitude of the gravitational acceleration, and $e z$ is a unit vector pointing vertically upward.
FIG. 1.

Schematic diagram of a droplet of fluid 2 falling in an unbounded oscillating fluid 1.

FIG. 1.

Schematic diagram of a droplet of fluid 2 falling in an unbounded oscillating fluid 1.

Close modal
Here, X and Y are in the horizontal plane, while Z is pointing vertically upward. The driving forcing on the droplet can be of different kinds, but one of the simplest, yet quite informative to consider, is one in which the domain is being periodically oscillated in the vertical direction, such that the position of the domain is moving vertically upward by a distance $A c H ( T )$ compared to the stationary reference frame. Here, Ac is the constant amplitude of the oscillation, while H(T) is the temporal part of the motion with a maximum value of unity. Hence, the far-field boundary condition is
$U → A c d H d T e z as R → ∞,$
(3)
where $A c d H / d T$ is the time-dependent forcing velocity and $R = X 2 + Y 2 + Z 2$ is the radial distance from the origin.
We let $F ̃ ( X , Y , Z , T ) = 0$ denote the equation of the interface between the two fluids. We require that the velocities are continuous at the interface, i.e.,
$U = U ̂ on F ̃ = 0.$
(4)
The kinematic condition, which implies that particles on the interface will remain on the interface, can be obtained by noting that since $F ̃$ is a scalar function which always vanishes at any point on the interface, its time derivative following any material point there should also vanish, i.e., $D F ̃ / D T = 0$, where D/DT is the convective derivative, or
$∂ F ̃ ∂ T + U · ∇ F ̃ = 0 on F ̃ = 0.$
(5)
The stress balance equation on the interface is
$( P ̂ − P ) n + ( τ ̃ ¯ − τ ̃ ̂ ¯ ) n = γ n ( ∇ · n ) − ∇ γ on F ̃ = 0,$
where $τ ̃ ¯ = 2 μ e ̃ ¯$ is the deviatoric stress tensor, with $e ̃ ¯ = 1 2 ( ∇ U + ( ∇ U ) T )$ the rate-of-strain tensor, γ is the surface tension, n is the unit outward pointing normal vector to the interface, and $μ = ν ρ$ is the dynamic viscosity. We shall assume that the surface tension is constant, and so, the tangential and normal stress balances are
$( ( τ ̃ ¯ − τ ̃ ̂ ¯ ) n ) · t i = 0 on F ̃ = 0 ,$
(6)
$P ̂ − P + ( ( τ ̃ ¯ − τ ̃ ̂ ¯ ) n ) · n = γ ∇ · n on F ̃ = 0 ,$
(7)
where $t i$ is a tangential vector.
We shall now introduce a coordinate system moving with the droplet and non-dimensionalize our governing equations as follows:
$T = R c U c t , ( X , Y ) = R c ( x , y ) , Z = R c ( z − t + A c R c h ) , ( U , V , U ̂ , V ̂ ) = U c ( u , v , u ̂ , v ̂ ) , W = U c ( w − 1 + A c R c d h d t ) , W ̂ = U c ( w ̂ − 1 + A c R c d h d t ) , τ ¯ ̃ = μ U c R c τ ¯ , e ¯ ̃ = U c R c e ¯ , P = μ U c R c p − ρ g R c z − ρ U c 2 A c R c d 2 h d t 2 z , P ̂ = μ U c R c p ̂ − ρ ̂ g R c z − ρ ̂ U c 2 A c R c d 2 h d t 2 z ,$
where lower case letters without tildes are used to denote the dimensionless parameters, Rc is the average radius of the droplet, and Uc is the average speed of the droplet which we assume to be non-zero. The dimensional timescale is $R c / U c$. To keep the problem as general as possible we are not specified the forcing function however, we are limiting the forcing frequency to a lower bound of $O ( 2 π U c / R c )$ so that the analysis is valid. To avoid confusion, we write $H ( T ) = h ( t )$ as H is a function of the dimensional time T, while h is a function of the dimensionless time t. We note that the hydrostatic force and the oscillations in the vertical direction have been included in the pressure. The dimensionless version of the Navier–Stokes and continuity equations (1) and (2) is
$Re ( ∂ u ∂ t + ( u · ∇ ) u ) = − ∇ p + ∇ 2 u , ∇ · u = 0 ,$
(8)
$κ Re ( ∂ u ̂ ∂ t + ( u ̂ · ∇ ) u ̂ ) = − ∇ p ̂ + λ ∇ 2 u ̂ , ∇ · u ̂ = 0$
(9)
with
$Re = R c U c ν , κ = ρ ̂ ρ and λ = μ ̂ μ,$
where $Re$ is the Reynolds number. In this study, we shall assume that the Reynolds number is sufficiently small that it can be neglected, so that the Navier–Stokes and continuity equations (8) and (9) reduce down to the Stokes flow equations
$∇ p = ∇ 2 u , ∇ · u = 0 , ∇ p ̂ = λ ∇ 2 u ̂ , ∇ · u ̂ = 0.$
(10)
We have assumed that κ is not too large and λ is not too small so that the terms in Eq. (9) are of a similar order to the corresponding terms in Eq. (8); otherwise, this would invalidate the use of Stokes flow. The far-field condition equation (3) is now written as
$u → e z as r → ∞ ,$
(11)
where $r = x 2 + y 2 + z 2$. Now, we use $F ( x , y , z , t ) = 0$ as the equation for the interface. The continuity of velocity condition equation (4) becomes
$u = u ̂ on F = 0,$
(12)
and the kinematic condition equation (5) becomes
$∂ F ∂ t + u · ∇ F = 0 on F = 0.$
(13)
Finally, the tangential stress condition equation (6) is
$( ( τ ¯ − τ ̂ ¯ ) n ) · t i = 0 on F = 0,$
(14)
and the normal stress condition equation (7) as
$p ̂ − p + [ 1 Fr 2 + A c R c d 2 h d t 2 ] Re ( 1 − κ ) z + ( ( τ ¯ − τ ̂ ¯ ) n ) · n = 1 Ca ∇ · n on F = 0 ,$
(15)
where
$Ca = μ U c γ and Fr = U c g R c ,$
where $Ca$ is the capillary number and $Fr$ is the Froude number.
To solve the Stokes flow equations in Eq. (10), we shall use axisymmetric spherical polar coordinates and write the velocities as
$u = u r ( r , θ ) e r + u θ ( r , θ ) e θ , u ̂ = u ̂ r ( r , θ ) e r + u ̂ θ ( r , θ ) e θ ,$
where $e r$ and $e θ$ are the unit vectors in the r and θ directions, respectively. We recall that r is the distance measured from the origin and θ is the angle measured anticlockwise from the positive z axis. For the functions $u ̂ r ( r , θ )$ and $u ̂ θ ( r , θ )$ we require,
$u ̂ r and u ̂ θ are bounded at r = 0.$
(16)
Using $e z = cos ( θ ) e r − sin ( θ ) e θ$, the far-field condition (11) becomes
$u r → cos ( θ ) and u θ → − sin ( θ ) as r → ∞,$
(17)
and the continuity of velocity condition (12) becomes
$u r = u ̂ r and u θ = u ̂ θ on F = 0.$
(18)
We now suppose that
$F = r − 1 − ε f ̃ ( θ , t ),$
where ε is a small constant representing the amplitude of the deviation of the droplet from a spherical droplet and $f ̃$ is an unknown function to be determined. Then, F = 0, the equation for the droplet interface, is given by
$r = 1 + ε f ̃,$
(19)
and the kinematic condition (13) is written as
$ε f ̃ t − u r + 1 r ε f ̃ θ u θ = 0 on r = 1 + ε f ̃ .$
(20)
Now, the unit normal and tangential vectors to the interface, $F = 0$, are
$n = 1 N ∇ F = 1 N ( e r − ε f ̃ θ r e θ ) and t = n × e ϕ = − 1 N ( ε f ̃ θ r e r + e θ ) ,$
where $N = | ∇ F | = 1 + ε 2 r − 2 f ̃ θ 2$. Finally, the tangential stress condition (14) is written as
$ε r f ̃ θ ( e r r − e θ θ − λ e ̂ r r + λ e ̂ θ θ ) = ( r 2 − ε 2 f ̃ θ 2 ) ( λ e ̂ r θ − e r θ ) on r = 1 + ε f ̃ ,$
(21)
where the components of the rate-of-strain tensor are as follows:
$e r r = ∂ u r ∂ r , e θ θ = 1 r ∂ u θ ∂ θ + u r r , e r θ = r 2 ∂ ∂ r ( u θ r ) + 1 2 r ∂ u r ∂ θ .$
The normal stress condition (15) is given by
$1 r 2 Ca ( ∂ ∂ r ( r 2 N ) − ε sin ( θ ) ∂ ∂ θ ( f ̃ θ sin ( θ ) N ) ) + p − p ̂ = ( 1 Fr 2 + A c R c d 2 h d t 2 ) Re ( 1 − κ ) r cos ( θ ) + 2 N 2 ( e r r − λ e ̂ r r − 2 ε f ̃ θ r ( e r θ − λ e ̂ r θ ) ) + 2 ε 2 f ̃ θ 2 r 2 N 2 ( e θ θ − λ e ̂ r θ ) on r = 1 + ε f ̃ .$
(22)
By introducing $η = cos ( θ )$, we can express the velocity components in terms of the streamfunctions ψ and $ψ ̂$ for the flow fields outside and inside of the droplet, respectively, as
$u r = − ψ η r 2 , u θ = − ψ r r 1 − η 2 , u ̂ r = − ψ ̂ η r 2 , u ̂ θ = − ψ ̂ r r 1 − η 2 .$
The Stokes flow equations (10) are satisfied when the streamfunctions satisfy equation (A3), namely,
$E 4 ψ = E 4 ψ ̂ = 0,$
(23)
where
$E 2 ψ = ∂ 2 ψ ∂ r 2 + 1 − η 2 r 2 ∂ 2 ψ ∂ η 2,$
(24)
which is explained in the book by Leal.39 The boundedness condition at r = 0 for the velocity components, Eq. (16), means that
$ψ ̂ r 2 is bounded at r = 0.$
(25)
The far-field condition (17) becomes
$ψ η → − r 2 η and ψ r → r ( 1 − η 2 ) as r → ∞ .$
By integrating these expressions and using $Q 1 = 1 2 ( η 2 − 1 )$, we obtain
$ψ → − r 2 Q 1 as r → ∞ .$
(26)
We now let $f ̃ ( θ , t ) = f ( η , t )$, which makes the equation for the interface
$r = ( 1 + ε f ) e r$
(27)
in vector form. The continuity of velocity condition (18) becomes
$ψ η = ψ ̂ η and ψ r = ψ ̂ r on r = 1 + ε f .$
Expanding each term using a Taylor series about ε = 0 and linearizing in ε, we obtain
$ψ η + ε f ψ η r = ψ ̂ η + ε f ψ ̂ η r on r = 1 ,$
(28)
$ψ r + ε f ψ r r = ψ ̂ r + ε f ψ ̂ r r on r = 1.$
(29)
The kinematic condition (20) now becomes
$ε f t + ψ η r 2 + ε f η ψ r r 2 = 0 on r = 1 + ε f .$
Expanding about ε = 0 and linearizing in ε yield
$ε f t + ψ η + ε f ψ r η − 2 ε f ψ η + ε f η ψ r = 0 on r = 1.$
(30)
The tangential stress condition (21) is also linearized in ε, and we obtain
$e r θ − λ e ̂ r θ = [ e r r − e θ θ + λ ( e ̂ θ θ − e ̂ r r ) ] ε r 1 − η 2 f η$
on $r = 1 + ε f$. We can express the components of the rate-of-strain tensor in terms of the stream function as
$e r r = 2 r 3 ψ η − ψ r η r 2 , e θ θ = ψ r η r 2 + η ψ r r 2 ( 1 − η 2 ) − ψ η r 3 , e r θ = ψ r r 2 1 − η 2 − ψ r r 2 r 1 − η 2 + 1 − η 2 2 r 3 ψ η η .$
The tangential stress condition is written as
$ψ r r 2 − ψ r r 2 r − Q 1 r 3 ψ η η + ε f η ( Q 1 ′ r 3 ψ r + 6 Q 1 r 4 ψ η − 4 Q 1 r 3 ψ r η ) = λ [ ψ ̂ r r 2 − ψ ̂ r r 2 r − Q 1 r 3 ψ ̂ η η + ε f η ( Q 1 ′ r 3 ψ ̂ r + 6 Q 1 r 4 ψ ̂ η − 4 Q 1 r 3 ψ ̂ r η ) ]$
on $r = 1 + ε f$. Upon expanding about ε = 0 and linearizing in ε, we get
$λ [ 2 ψ ̂ r − ψ ̂ r r − 2 Q 1 ψ ̂ η η + ε f ( 3 ψ ̂ r r − 4 ψ ̂ r − ψ ̂ rrr ) ] + 2 λ ε f η ( Q 1 ′ ψ ̂ r + 6 Q 1 ψ ̂ η − 4 Q 1 ψ ̂ r η ) − 2 λ ε f Q 1 ( ψ ̂ r η η − 3 ψ ̂ η η ) = 2 ψ r − ψ r r − 2 Q 1 ψ η η + 2 ε f η ( Q 1 ′ ψ r + 6 Q 1 ψ η − 4 Q 1 ψ r η ) + ε f [ 3 ψ r r − 4 ψ r − ψ rrr − 2 Q 1 ( ψ r η η − 3 ψ η η ) ]$
(31)
on r = 1. We now turn to the normal stress condition (22), which is also linearized in ε to give
$2 ( e r r + 2 ε r 1 − η 2 f η e r θ ) − 2 λ ( e ̂ r r + 2 ε r 1 − η 2 f η e ̂ r θ ) + p ̂ − p + ( 1 Fr 2 + A c R c d 2 h d t 2 ) Re ( 1 − κ ) r η = 2 Ca ( 1 r + ε r 2 ( Q 1 f η ) η ) on r = 1 + ε f ,$
which upon substituting in the components of the rate-of-strain tensor yields
$2 Ca ( 1 r + ε r 2 ( Q 1 f η ) η ) = p ̂ − p + 4 r 3 ψ η − 2 ψ r η r 2 + 2 ε f η ( 2 r 3 ψ r − 1 r 2 ψ r r − 2 r 4 Q 1 ψ η η ) − 2 λ [ 2 r 3 ψ ̂ η − ψ ̂ r η r 2 + ε f η ( 2 r 3 ψ ̂ r − 1 r 2 ψ ̂ r r − 2 r 4 Q 1 ψ ̂ η η ) ] + ( 1 Fr 2 + A c R c d 2 h d t 2 ) Re ( 1 − κ ) r Q 1 ′$
on $r = 1 + ε f$.
To carry out this expansion, we are going to assume that $A c / R c$ and $Fr$ are large so that the leading terms are retained in the expansion. In particular, we choose $A c / R c = O ( ε − 1 ) , Fr = O ( ε − 1 )$, and $Re = O ( ε 2 )$, with f and h both of O(1). Then, expanding about ε = 0 and linearizing, we obtain
$2 Ca ( 1 − ε f + ε ( Q 1 f η ) η ) = p ̂ − p + ε f ( p ̂ r − p r ) + ( 1 + ε f Fr 2 + A c R c d 2 h d t 2 ) Re ( 1 − κ ) Q 1 ′ + 4 ψ η − 2 ψ r η + 2 ε [ f ( 4 ψ r η − 6 ψ η − ψ r r η ) + f η ( 2 ψ r − ψ r r − 2 Q 1 ψ η η ) ] − 2 λ ε f ( 4 ψ ̂ r η − 6 ψ ̂ η − ψ ̂ r r η ) − 2 λ [ 2 ψ ̂ η − ψ ̂ r η ] − 2 λ ε f η ( 2 ψ ̂ r − ψ ̂ r r − 2 Q 1 ψ ̂ η η ) on r = 1 .$
(32)
To solve the above system of equations, we introduce the expansions
$ψ ( r , η , t ) = ψ 0 ( r , η ) + ε ψ 1 ( r , η , t ) , ψ ̂ ( r , η , t ) = ψ ̂ 0 ( r , η ) + ε ψ ̂ 1 ( r , η , t ) , p ( r , η , t ) = p 0 ( r , η ) + ε p 1 ( r , η , t ) , p ̂ ( r , η , t ) = p ̂ 0 ( r , η ) + ε p ̂ 1 ( r , η , t ) .$
We begin by seeking the zeroth-order solutions, which are independent of time and the droplet is a sphere of radius 1. The zeroth-order Stokes flow equations are
$E 4 ψ 0 = 0 and E 4 ψ ̂ 0 = 0,$
whose general solutions are given by Eqs. (B3) and (B4), namely,
$ψ 0 = ∑ j = 1 ∞ ( A j 0 r j + 3 + C j 0 r j + 1 + B j 0 r 2 − j + D j 0 r − j ) Q j , ψ ̂ 0 = ∑ j = 1 ∞ ( A ̂ j 0 r j + 3 + C ̂ j 0 r j + 1 + B ̂ j 0 r 2 − j + D ̂ j 0 r − j ) Q j ,$
where the Qjs are a modified set of Gegenbauer polynomials, which satisfy Eq. (C1) in  Appendix C. Additional properties of the modified Gegenbauer polynomials are given in  Appendix C. The boundedness condition for the velocity components r = 0 (25) yields
$ψ ̂ 0 r 2 is bounded at r = 0.$
This means that $B ̂ j 0 = D ̂ j 0 = 0$ for $j ≥ 1$. The far-field condition (26) then yields
$ψ 0 → − r 2 Q 1 as r → ∞ .$
This means that $A j 0 = 0$ for $j ≥ 1 , C 1 0 = − 1$, and $C j 0 = 0$ for $j ≥ 2$. The kinematic condition (30) yields
$ψ η 0 = 0 on r = 1.$
As the derivatives of the Gegenbauer polynomials are linearly independent, we can equate their coefficients to yield $D 1 0 = 1 − B 1 0$ and $D j 0 = − B j 0$ for $j ≥ 2$. Next, the continuity of velocity conditions (28) and (29) gives
$ψ η 0 = ψ ̂ η 0 and ψ r 0 = ψ ̂ r 0 on r = 1.$
The first condition gives $A ̂ j 0 = − C ̂ j 0$ for $j ≥ 1$. The second gives $A ̂ 1 0 = − β / 2$, where $β = 3 − 2 B 1 0$ and $A ̂ j 0 = B j 0$ for $j ≥ 2$. We next turn to the tangential stress condition (31), which is written as
$λ [ 2 ψ ̂ r 0 − ψ ̂ r r 0 − 2 Q 1 ψ ̂ η η 0 ] = 2 ψ r 0 − ψ r r 0 − 2 Q 1 ψ η η 0$
on r = 1. This equation gives
$B 1 0 = 2 + 3 λ 2 ( 1 + λ )$
(33)
and $B j 0 = 0$ for $j ≥ 2$. At this stage, the leading-order stream functions are given by
$ψ 0 = ( B 1 0 r − r 2 + 1 − B 1 0 r ) Q 1 ,$
(34)
$ψ ̂ 0 = β 2 ( r 2 − r 4 ) Q 1 .$
(35)
Since $η = cos ( θ )$ and $Q 1 = 1 2 ( η 2 − 1 )$, we can write
$ψ 0 = ( r 2 2 − B 1 0 r 2 + B 1 0 − 1 2 r ) sin 2 ( θ ) , ψ ̂ 0 = β 4 ( r 4 − r 2 ) sin 2 ( θ ) .$
As λ is a positive constant, $1 ≤ B 1 0 ≤ 3 2$, and so, $ψ 0$ is a non-negative monotonically increasing function of r, while $ψ ̂ 0$ is a non-positive function of r, which has a local minimum at $r = 1 / 2$, with the minimum value $ψ ̂ 0 = − β / 16$. This means that a recirculation zone exists inside the droplet. We notice that if $λ → ∞$, then $B 1 0 → 3 / 2$ and $ψ ̂ 0 → 0$. Figure 2 illustrates the streamlines obtained using Eqs. (34) and (35). The general solutions for $p 0$ and $p ̂ 0$ are given by Eqs. (B5) and (B6), namely,
$p 0 = p 0 + ∑ j = 1 ∞ ( − A j 0 4 j + 6 j r j + B j 0 2 − 4 j j + 1 r − j − 1 ) Q j ′$
(36)
and
$p ̂ 0 = p ̂ 0 + λ ∑ j = 1 ∞ ( − A ̂ j 0 4 j + 6 j r j + B ̂ j 0 2 − 4 j j + 1 r − j − 1 ) Q j ′$
(37)
as presented in  Appendix B and shown by in the book Leal.39 Substituting in the values obtained for the coefficients $A j 0 , B j 0 , A ̂ j 0$, and $B ̂ j 0$ into Eqs. (36) and (37) yields
$p 0 = p 0 − B 1 0 r 2 Q 1 ′ and p ̂ 0 = p ̂ 0 + 5 λ β r Q 1 ′ .$
(38)
The normal stress condition (32) yields
$2 Ca = p ̂ 0 − p 0 + Re Fr 2 ( 1 − κ ) Q 1 ′ − λ [ 4 ψ ̂ η 0 − 2 ψ ̂ r η 0 ] + 4 ψ η 0 − 2 ψ r η 0 on r = 1 .$
(39)
Substituting in the solutions for $p ̂ 0 , p 0$, $ψ 0$, and $ψ ̂ 0$ into Eq. (39) and equating coefficients of the constants and $Q 1 ′$ yields
$2 Ca = p ̂ 0 − p 0 and 3 B 1 0 = Re Fr 2 ( κ − 1 ) .$
(40)
Finally, using the definitions of the dimensionless parameters, the second equation above becomes
$6 + 9 λ 2 ( 1 + λ ) = ρ g R c 2 μ U c ( κ − 1 ) .$
Substituting in the ratios for κ and λ yields
$U c = 2 g R c 2 ( ρ ̂ − ρ ) ( μ + μ ̂ ) 3 μ ( 2 μ + 3 μ ̂ ) .$
(41)
Hence, the speed of the droplet is proportional to gravity, the difference in density, and the square of the droplet radius, which are the same results obtained by Leal.39
FIG. 2.

Contours of the leading-order streamfunctions $ψ 0$ and $ψ ̂ 0$ in Eqs. (34) and (35) with $λ = 1.3$. The contour values for $ψ 0$ are $7 × 10 − 5$, 0.07, 0.35, 1, 2 and $ψ ̂ 0$ are $− 1 × 10 − 6$, −0.000 025, −0.003, −0.015, and −0.026. The red and blue lines represent the streamlines for $ψ ̂ 0$ and $ψ 0$, respectively.

FIG. 2.

Contours of the leading-order streamfunctions $ψ 0$ and $ψ ̂ 0$ in Eqs. (34) and (35) with $λ = 1.3$. The contour values for $ψ 0$ are $7 × 10 − 5$, 0.07, 0.35, 1, 2 and $ψ ̂ 0$ are $− 1 × 10 − 6$, −0.000 025, −0.003, −0.015, and −0.026. The red and blue lines represent the streamlines for $ψ ̂ 0$ and $ψ 0$, respectively.

Close modal
We next turn our attention to the first-order solutions in which the droplet interface is allowed to deform in η and time. The first-order Stokes flow equations are
$E 4 ψ 1 = 0 and E 4 ψ ̂ 1 = 0,$
whose general solutions are
$ψ 1 = ∑ j = 1 ∞ ( A j 1 r j + 3 + C j 1 r j + 1 + B j 1 r 2 − j + D j 1 r − j ) Q j , ψ ̂ 1 = ∑ j = 1 ∞ ( A ̂ j 1 r j + 3 + C ̂ j 1 r j + 1 + B ̂ j 1 r 2 − j + D ̂ j 1 r − j ) Q j .$
The boundedness condition for the velocities at r = 0 (25) yields
$ψ ̂ 1 r 2 is bounded at r = 0 ,$
which means that $B ̂ j 1 = D ̂ j 1 = 0$ for $j ≥ 1$. The far-field condition (26) yields
$ψ 1 r 2 → 0 as r → ∞,$
which means that $A j 1 = C j 1 = 0$ for $j ≥ 1$. Turning then to the continuity of velocity conditions (28) and (29),
$ψ η 1 + f ψ η r 0 = ψ ̂ η 1 + f ψ ̂ η r 0 on r = 1 , ψ r 1 + f ψ r r 0 = ψ ̂ r 1 + f ψ ̂ r r 0 on r = 1.$
Substituting in $ψ 0$ and $ψ ̂ 0$ reduces these to
$ψ η 1 = ψ ̂ η 1 and ψ r 1 = ψ ̂ r 1 + α f Q 1 on r = 1 ,$
where $α = 12 B 1 0 − 15$. As the Qjs are linearly independent functions, the first condition implies that $B j 1 = A ̂ j 1 + C ̂ j 1 − D j 1$ for $j ≥ 1$. Using this result, the second condition gives
$f = 1 Q 1 ∑ j = 1 ∞ Λ j Q j,$
(42)
where
$Λ j = − 1 α ( ( 2 j + 1 ) A ̂ j 1 + ( 2 j − 1 ) C ̂ j 1 + 2 D j 1 ),$
(43)
for $j ≥ 1$. Here, we are assuming that $α ≠ 0$, i.e., $B 1 0 ≠ 5 4$, i.e., $λ ≠ 1$, i.e., $μ ≠ μ ̂$ and the two fluids have different dynamic viscosities. We notice that at this stage, we can write the first-order streamfunctions as
$ψ 1 = ∑ j = 1 ∞ ( ( A ̂ j 1 + C ̂ j 1 − D j 1 ) r 2 − j + D j 1 r − j ) Q j ,$
(44)
$ψ ̂ 1 = ∑ j = 1 ∞ ( A ̂ j 1 r j + 3 + C ̂ j 1 r j + 1 ) Q j .$
(45)
Turning now to the kinematic condition (30), we have
$f t + ψ η 1 + f ψ r η 0 − 2 f ψ η 0 + f η ψ r 0 = 0 on r = 1.$
Substituting in $ψ 0$ reduces this to
$f t + ψ η 1 − β ( Q 1 f ) η = 0 on r = 1 ,$
from which
$∑ j = 1 ∞ ( d Λ j d t Q j + Ω j Q j ′ Q 1 ) = 0,$
where
$Ω j = A ̂ j 1 + C ̂ j 1 − β Λ j = A ̂ j 1 + C ̂ j 1 + β α ( ( 2 j + 1 ) A ̂ j 1 + ( 2 j − 1 ) C ̂ j 1 + 2 D j 1 ) ,$
(46)
for $j ≥ 1$. We now make use of the identity (C5) in  Appendix C to obtain
$∑ j = 1 ∞ d Λ j d t Q j = ∑ j = 1 ∞ Ω j + 1 j ( j + 1 ) Q j 4 j + 6 − ∑ j = 2 ∞ Ω j − 1 j ( j + 1 ) Q j 4 j − 2 .$
By defining $Ω 0 = 0$ and since the Qjs are linearly independent, we obtain
$d Λ j d t = j ( j + 1 ) ( Ω j + 1 4 j + 6 − Ω j − 1 4 j − 2 ) for j ≥ 1.$
(47)
We now turn to the tangential stress condition (31), and substituting in $ψ 0$ and $ψ ̂ 0$ leads to
$λ [ 2 ψ ̂ r 1 − ψ ̂ r r 1 − 2 Q 1 ψ ̂ η η 1 ] + 2 α Q 1 Q 1 ′ f η = 2 ψ r 1 − ψ r r 1 − 2 Q 1 ψ η η 1 + 15 ( 1 − B 1 0 ) Q 1 f on r = 1 ,$
which becomes
$∑ j = 1 ∞ ( j 2 ( 1 − λ ) − ( 1 + 3 λ ) j − 2 ) A ̂ j 1 Q j ∑ j = 1 ∞ [ ( 1 − λ ) ( j 2 − j − 2 ) C ̂ j 1 + ( 4 j + 2 ) D j 1 ] Q j + 2 ( 1 − λ ) ∑ j = 1 ∞ ( A ̂ j 1 + C ̂ j 1 ) Q 1 Q j ″ = − 2 α Q 1 Q 1 ′ f η + 15 ( 1 − B 1 0 ) Q 1 f .$
Using Eq. (C1), we obtain
$∑ j = 1 ∞ Ψ j Q j = α Q 1 Q 1 ′ f η + 15 2 ( B 1 0 − 1 ) Q 1 f ,$
where we have defined
$Ψ j = ( j 2 ( λ − 1 ) + 2 λ j + 1 ) A ̂ j 1 + ( λ − 1 ) ( j 2 − 1 ) C ̂ j 1 − ( 2 j + 1 ) D j 1 for j ≥ 1 .$
(48)
Next, we use the shape of the droplet, Eq. (42), along with the identity (C7) in  Appendix C to get
$∑ j = 1 ∞ Γ j Q j = α ∑ j = 1 ∞ Λ j ( ∑ m = 1 j − 2 ( − 1 ) j − m + 1 2 j ( j + 1 ) m ( m + 1 ) ( 2 m + 1 ) Q m ) ,$
where we have defined as
$Γ j = Ψ j − ( α j + 1 2 ( 15 − 9 B 1 0 ) ) Λ j = ( j 2 ( λ + 1 ) + ( 2 λ + 1 ) j + 1 ) A ̂ j 1 + ( j 2 ( λ + 1 ) − j + 1 − λ ) C ̂ j 1 − D j 1 + 15 − 9 B 1 0 2 α ( ( 2 j + 1 ) A ̂ j 1 + ( 2 j − 1 ) C ̂ j 1 + 2 D j 1 ) ,$
(49)
for $j ≥ 1$. By changing the order of the summations and using the fact that the Qjs are linearly independent functions, the tangential stress condition reduces to
$Γ j = α Φ j for j ≥ 1,$
(50)
where
$Φ j = j ( j + 1 ) ( 2 j + 1 ) ∑ m = j + 2 ∞ Λ m ( ( − 1 ) m − j + 1 ) 2 m ( m + 1 ),$
for $j ≥ 1$. By setting $m = 2 i + j$, the above expression is simplified to
$Φ j = ∑ i = 1 ∞ j ( j + 1 ) ( 2 j + 1 ) ( 2 i + j ) ( 2 i + j + 1 ) Λ 2 i + j for j ≥ 1.$
(51)
We next look at the normal stress condition (32) and substituting in $ψ 0$ and $ψ ̂ 0$, from Eqs. (34) and (35) and p0 and $p ̂ 0$ from Eq. (38), and using Eq. (40) lead to
$2 Ca ( ( Q 1 f η ) η − f ) = p ̂ 1 − p 1 − ξ Q 1 ′ f − 3 B 1 0 d ζ d t Q 1 ′ − λ [ 4 ψ ̂ η 1 − 2 ψ ̂ r η 1 ] + 4 ψ η 1 − 2 ψ r η 1 on r = 1 ,$
where $ξ = 30 − 21 B 1 0$ and $ζ = A c Fr 2 ε R c d h d t$. Using the general solutions for the pressure, we have
$p 1 = ∑ j = 1 ∞ ( A ̂ j 1 + C ̂ j 1 − D j 1 ) ( 2 − 4 j ) Q j ′ ( j + 1 ) r j + 1 ,$
(52)
$p ̂ 1 = − λ ∑ j = 1 ∞ A ̂ j 1 4 j + 6 j r j Q j ′ .$
(53)
Substituting the pressures into our normal stress condition yields
$2 Ca ( ( Q 1 f η ) η − f ) = ∑ j = 1 ∞ 2 Θ j Q j ′ − ξ Q 1 ′ f − 3 B 1 0 d ζ d t Q 1 ′,$
where we have defined as
$Θ j = ( λ ( j − 1 ) + j + 2 j − 1 j + 1 ) ( A ̂ j 1 + C ̂ j 1 ) − 3 λ j A ̂ j 1 + 3 D j 1 j + 1 for j ≥ 1.$
(54)
We now use the shape of the droplet in Eq. (42) to get
$2 Ca ∑ j = 1 ∞ Λ j ( Q 1 ′ Q 1 ( Q j Q 1 ) ′ + Q 1 2 ( Q j Q 1 ) ″ − Q j ) = ∑ j = 1 ∞ 2 Θ j Q j ′ Q 1 − ξ Q 1 ′ ∑ j = 1 ∞ Λ j Q j − 3 B 1 0 d ζ d t Q 1 ′ Q 1 .$
To simplify this expression, we make use of the identities (C4)–(C7) in  Appendix C to give
$2 Ca ∑ j = 1 ∞ Λ j ( 1 2 ( j − 2 ) ( j + 1 ) Q j ) − 2 Ca ∑ j = 1 ∞ Λ j ( ∑ m = 1 j − 2 ( − 1 ) j − m + 1 2 j ( j + 1 ) m ( m + 1 ) ( 2 m + 1 ) Q m ) = ∑ j = 1 ∞ Θ j 2 j + 1 ( ( j + 1 ) ( j + 2 ) Q j + 1 − j ( j − 1 ) Q j − 1 ) − ξ ∑ j = 1 ∞ Λ j 2 j + 1 ( ( j + 2 ) Q j + 1 + ( j − 1 ) Q j − 1 ) − 3 B 1 0 d ζ d t Q 2 .$
As before, we change the order of the summations and simplify the above equation to
$1 Ca ∑ j = 1 ∞ ( ( j − 2 ) ( j + 1 ) Λ j − 2 Φ j ) Q j = − 3 B 1 0 d ζ d t Q 2 + ∑ j = 1 ∞ j ( j + 1 ) Q j ( Θ j − 1 2 j − 1 − Θ j + 1 2 j + 3 ) − ξ ∑ j = 1 ∞ ( ( j + 1 ) Λ j − 1 2 j − 1 + j Λ j + 1 2 j + 3 ) Q j ,$
where we have defined $Λ 0 = Θ 0 = 0$. As the Qjs are linearly independent, the normal stress condition reduces to
$1 Ca ( ( j − 2 ) ( j + 1 ) Λ j − 2 Γ j α ) = − 3 B 1 0 d ζ d t δ 2 , j + j ( j + 1 ) ( Θ j − 1 2 j − 1 − Θ j + 1 2 j + 3 ) − ξ ( ( j + 1 ) Λ j − 1 2 j − 1 + j Λ j + 1 2 j + 3 ) for j ≥ 1 ,$
(55)
where we have substituted (50) in the tangential stress condition, and $δ i , j$ is the Kronecker delta, equal to 1 when i = j and zero otherwise.
Equations (43), (46), (49), and (54) are expressions for Λj, Ωj, $Γ j$, and Θj in terms of $A ̂ j 1 , C ̂ j 1$, and $D j 1$. We notice that we can invert Eqs. (43), (46), and (49) to yield
$A ̂ j 1 = 1 2 ( 6 B 1 0 − 7 − 2 j 2 ) a j + α β 2 j ( j + 1 ) c j − β j ( j + 1 ) ( 2 j 2 + 15 B 1 0 − 23 ) 4 ( 2 j + 1 ) b j ,$
(56)
$C ̂ j 1 = 1 2 ( 2 j 2 + 4 j + 3 β ) a j − α β 2 j ( j + 1 ) c j + β j ( j + 1 ) ( 2 j 2 + 4 j − 21 + 15 B 1 0 ) 4 ( 2 j + 1 ) b j ,$
(57)
$D j 1 = − j ( j + 1 ) ( 2 β j 2 + 2 α j + 30 ( B 1 0 ) 2 − 77 B 1 0 + 51 ) 4 ( 2 j + 1 ) b j + ( 4 − 3 B 1 0 − j 2 ) a j − α β 2 j ( j + 1 ) c j .$
(58)
To express the system of Eqs. (47), (50), (51), and (55) in a simpler form, we define the new variables
$a j = Ω j 2 j + 1 , b j = 2 Λ j j ( j + 1 ) ,$
(59)
$c j = 2 Γ j α j ( j + 1 ) ( 2 j + 1 ) and d j = 2 Θ j 2 j + 1,$
(60)
where $a 0 = b 0 = d 0 = 0$. Furthermore, using the new variables, we can write
$d j = q j a a j + q j b b j + q j c c j,$
(61)
where
$q j a = 6 ( 2 j + 3 B 1 0 − 2 ) j ( j + 1 ) ( 2 j + 1 ) + 4 j 3 + 4 j 2 − 2 j − 6 β j ( 2 j + 1 ) ,$
(62)
$q j b = 4 j 4 + 8 j 3 + ( 92 − 72 B 1 0 ) j 2 2 ( 2 j + 1 ) 2 + ( 3 B 1 0 − 17 ) j + 6 ( B 1 0 − 1 ) ( 15 B 1 0 − 23 ) 2 ( 2 j + 1 ) 2 ,$
(63)
$q j c = − 3 α ( j + 2 B 1 0 − 2 ) 2 j + 1,$
(64)
for $j ≥ 1$ with $q 0 a = q 0 b = q 0 c = 0$. The kinematic condition equation (47) becomes
$d b j d t = a j + 1 − a j − 1 for j ≥ 1.$
(65)
Substituting Eq. (51) into the tangential stress condition equation (50) with these new variables yields
$c j = ∑ i = 1 ∞ b 2 i + j for j ≥ 1.$
(66)
Next, the normal stress condition equation (55) becomes
$1 Ca ( ( j − 2 ) ( j + 1 ) b j − c j ( 4 j + 2 ) ) = − B 1 0 d ζ d t δ 2 , j + d j − 1 − d j + 1 − ξ ( j − 1 2 j − 1 b j − 1 + j + 2 2 j + 3 b j + 1 )$
for $j ≥ 1$. We can eliminate dj from this equation to obtain
$1 Ca ( ( j − 2 ) ( j + 1 ) b j − c j ( 4 j + 2 ) ) = − B 1 0 d ζ d t δ 2 , j + q j − 1 a a j − 1 + q j − 1 b b j − 1 + q j − 1 c c j − 1 − q j + 1 a a j + 1 − q j + 1 b b j + 1 − q j + 1 c c j + 1 − ξ ( j − 1 2 j − 1 b j − 1 + j + 2 2 j + 3 b j + 1 ) for j ≥ 1.$
(67)
Hence, our system of equations consists now of Eqs. (65)–(67) which we need to solve to obtain aj, bj, and cj. Once bj is known, we can obtain Λj from Eq. (59), and then, using Eq. (42), we can obtain f, and then, from $r = 1 + ε f$, we can obtain the shape of the droplet. Using Eqs. (56)–(58), we can obtain $A ̂ j 1 , C ̂ j 1$, and $D j 1$. The streamfunctions and the pressures can then be obtained from Eqs. (44), (45), and (53).
The drag on the droplet only acts in the z direction and is given by
$D f = ∬ τ ¯ ̃ n · e z d S ̃ = μ U c R c C d,$
where the drag coefficient is given by
$C d = ∬ τ ¯ n · e z d S .$
(68)
In  Appendix D, we find that the drag coefficient is given by Eq. (68), namely,
$C d = 4 π B 1 0 + 4 3 ε π ( A ̂ 1 1 + C ̂ 1 1 − D 1 1 ) = 4 π B 1 0 + 4 3 ε π ( 3 B 1 0 a 1 + α β c 1 ) + 2 3 ε π ( 15 − 23 B 1 0 + 10 ( B 1 0 ) 2 ) b 1 .$
Notice, the relation for the drag coefficient reduces to the one for the unforced case that of a free droplet given by Leal.39 The volume of the droplet is given by
$V = ∫ 0 2 π ∫ 0 π ∫ 0 1 + ε f r 2 sin ( θ ) drd θ d ϕ = 2 π 3 ∫ 0 π ( 1 + ε f ) 3 sin ( θ ) d θ = 2 π 3 ∫ − 1 1 ( 1 + ε f ) 3 d η .$
Linearizing in ε gives
$V ≈ 2 π 3 ∫ − 1 1 1 + 3 ε f d η = 4 π 3 + 2 π ε ∫ − 1 1 f d η .$
Using Eqs. (42) and (C9), we obtain
$∫ − 1 1 f d η = ∫ − 1 1 1 Q 1 ∑ j = 1 ∞ Λ j Q j d η = ∑ j = 1 ∞ 2 Λ 2 j − 1 j ( 2 j − 1 ) = ∑ j = 1 ∞ 2 b 2 j − 1 .$
Conservation of mass dictates that for the droplet to have a volume of $4 π / 3$, corresponding to the volume of the undisturbed droplet, a sphere of radius 1, we require that $∫ − 1 1 f d η = 0$, or
$∑ j = 1 ∞ b 2 j − 1 = 0.$
(69)
We note that as the droplet is falling and that the domain is being forced to move in the vertical direction, the center of the mass of the droplet will also be moving vertically in the domain. We shall now consider how the center of mass changes in our coordinate system. The center of mass of the droplet is at $( 0 , 0 , z c )$, where zc is given by
$z c = ∫ 0 2 π ∫ 0 π ∫ 0 1 + ε f z r 2 sin ( θ ) drd θ d ϕ ∫ 0 2 π ∫ 0 π ∫ 0 1 + ε f r 2 sin ( θ ) drd θ d ϕ = 3 ∫ − 1 1 ( 1 + ε f ) 4 η d η 4 ∫ − 1 1 ( 1 + ε f ) 3 d η ≈ 3 ∫ − 1 1 ( 1 + 4 ε f ) η d η 4 ∫ − 1 1 1 + 3 ε f d η = 3 ε 2 ∫ − 1 1 f η d η = 3 ε 2 ∑ j = 1 ∞ Λ j ∫ − 1 1 Q 1 ′ Q 1 Q j d η$
using $∫ − 1 1 f d η = 0$, so that the droplet volume remains as $4 π / 3$. Using Eq. (C10), we obtain
$z c ≈ 3 ε ∑ k = 1 ∞ Λ 2 k k ( 2 k + 1 )$
(70)
for the vertical coordinate of the center of the droplet.
Before considering the case when the droplet is forced to oscillate, in this section, we consider the case when the droplet is falling and allowed to deform. For a steady flow field, Eq. (65) reduces to
$a j + 1 = a j − 1 for j ≥ 1.$
This means that $a 2 j = 0$ and $a 2 j − 1 = a 1$ for $j ≥ 1$. However, if the $a 2 j − 1$ s are all equal, then this means that the $Ω 2 j − 1$ s tend to infinity as j tends to infinity, which means that at least one of the $A ̂ j$ s, $C ̂ j$ s, or Dj s would diverge. Hence, in order to obtain a convergent solution, we require that $a 1 = 0$. Equation (66) reduces to
$c j = ∑ i = 1 ∞ b 2 i + j for j ≥ 1.$
Substituting this into Eq. (67) yields
$1 Ca ( ( j − 2 ) ( j + 1 ) b j − ( 4 j + 2 ) ∑ i = 1 ∞ b 2 i + j ) = q j − 1 b b j − 1 − q j + 1 b b j + 1 + q j − 1 c ∑ i = 1 ∞ b 2 i + j − 1 − q j + 1 c ∑ i = 1 ∞ b 2 i + j + 1 − ξ ( j − 1 2 j − 1 b j − 1 + j + 2 2 j + 3 b j + 1 ) for j ≥ 1.$
(71)
Using Eq. (69) for a steady state, we have
$∑ i = 1 ∞ b 2 i + 2 j − 1 = ∑ i = 1 ∞ b 2 i − 1 − ∑ i = 1 j b 2 i − 1 = − ∑ i = 1 j b 2 i − 1 .$
For convenience, we introduce the constant I0 defined as
$I 0 = ∑ i = 1 ∞ b 2 i so ∑ i = 1 ∞ b 2 i + 2 j = ∑ i = 1 ∞ b 2 i − ∑ i = 1 j b 2 i = I 0 − ∑ i = 1 j b 2 i .$
Then, using the odd values in Eq. (71) and simplifying we have
$2 Ca ( j ( 2 j − 3 ) b 2 j − 1 + ( 4 j − 1 ) ∑ i = 1 j b 2 i − 1 ) = ( q 2 j − 2 b − 2 j − 2 4 j − 3 ξ ) b 2 j − 2 + ( q 2 j − 2 c − q 2 j c ) I 0 − ( q 2 j b − q 2 j − 2 c + 2 j + 1 4 j + 1 ξ ) b 2 j + ( q 2 j c − q 2 j − 2 c ) ∑ i = 1 j b 2 i$
for $j ≥ 1$. Next, using the even values in Eq. (71) and simplifying we have
$2 Ca ( ( j − 1 ) ( 2 j + 1 ) b 2 j − ( 4 j + 1 ) ( I 0 − ∑ i = 1 j b 2 i ) ) = ( q 2 j − 1 b − 2 j − 1 4 j − 1 ξ ) b 2 j − 1 + ( q 2 j + 1 c − q 2 j − 1 c ) ∑ i = 1 j b 2 i − 1 + ( q 2 j + 1 c − q 2 j + 1 b − 2 j + 2 4 j + 3 ξ ) b 2 j + 1$
for $j ≥ 1$. To solve these in the small-capillary-number limit, we let $b 2 = χ$, where χ is an undetermined constant that measures the deformation of the surface of the droplet. We can now obtain a non-trivial solution and use the expansions
$I 0 = χ + Ca 2 I 0 ( 2 ) + O ( Ca 4 ) , b 2 j − 1 = Ca b 2 j − 1 ( 1 ) + O ( Ca 3 ) for j ≥ 1 , b 2 j = χ δ 1 , j + Ca 2 b 2 j ( 2 ) + O ( Ca 4 ) for j ≥ 1.$
Substituting these into the odd equations at leading order gives
$2 j ( 2 j − 3 ) b 2 j − 1 ( 1 ) + 2 ( 4 j − 1 ) ∑ i = 1 j b 2 i − 1 ( 1 ) = χ ( q 2 b − 2 5 ξ ) δ 2 , j − χ ( q 2 b + 3 5 ξ ) δ 1 , j for j ≥ 1.$
Solving this yields
$b 1 ( 1 ) = − χ 4 ( q 2 b + 3 5 ξ ) , b 3 ( 1 ) = χ 4 ( q 2 b + 17 45 ξ ) , b 2 j − 1 ( 1 ) = − ( 4 j − 1 ) ( 2 j − 1 ) ( j + 1 ) ∑ i = 1 j − 1 b 2 i − 1 ( 1 ) for j ≥ 2.$
(72)
Solving the recursion equation in Eq. (72) gives
$b 2 j − 1 ( 1 ) = ( 4 j − 1 ) χ ξ 2 j ( j + 1 ) ( 2 j − 1 ) ( 2 j − 3 ) for j ≥ 3.$
Using now Eq. (66), we have
$c 2 j = O ( Ca 2 ) and c 2 j − 1 = − ∑ i = 1 j b 2 i − 1 for j ≥ 1.$
Hence,
$c 1 = χ Ca 4 ( q 2 b + 3 5 ξ ) + O ( Ca 3 ) , c 3 = χ ξ 18 Ca + O ( Ca 3 ) , c 2 j − 1 = χ Ca ξ 2 ( j + 1 ) ( 2 j − 1 ) + O ( Ca 3 ) for j ≥ 3.$
Using Eq. (68) for the drag coefficient, we have
$C d = 4 π B 1 0 − ε χ π 10 Ca ( 3 ( B 1 0 ) 2 − 38 B 1 0 + 50 ) ( 58 ( B 1 0 ) 2 − 155 B 1 0 + 105 ) .$
We notice that when $χ > 0$, the $O ( Ca )$ term is reducing the drag coefficient as long as
$B 1 0 < 19 − 211 3 ≈ 1.4914 ,$
which corresponds to $λ = 28 + 2 211 ≈ 57.05$. Using Eq. (59), we obtain
$Λ 1 = − χ Ca 4 ( q 2 b + 3 5 ξ ) + O ( Ca 3 ) , Λ 2 = 3 χ , Λ 3 = 3 χ Ca 2 ( q 2 b + 17 45 ξ ) + O ( Ca 3 ) , Λ 2 j = O ( Ca 2 ) for j ≥ 2 , Λ 2 j − 1 = χ Ca ξ ( 4 j − 1 ) 2 ( j + 1 ) ( 2 j − 3 ) + O ( Ca 3 ) for j ≥ 3.$
Using Eq. (70), we find that the vertical coordinate of the center of mass of the droplet is given by $z c = 3 ε χ + O ( Ca 2 )$. Furthermore, from Eq. (42), we obtain
$f = 3 χ Q 2 Q 1 − χ Ca 4 ( q 2 b + 3 5 ξ ) + 3 χ Ca 2 ( q 2 b + 17 45 ξ ) Q 3 Q 1 + ∑ j = 3 ∞ χ Ca ξ ( 4 j − 1 ) 2 ( j + 1 ) ( 2 j − 3 ) Q 2 j − 1 Q 1 + O ( Ca 2 ) .$
(73)
Using $x = r cos ( θ )$ and $z = r sin ( θ )$ with $η = cos ( θ )$ and $r = 1 + ε f$, we can construct the droplet interface. For convenience, we truncated the series in Eq. (73) to the first 25 terms. Typical droplet interface profiles are displayed in Fig. 3 using $ε = 1 16$ and $Ca = 1 4$ for various values of λ with (a) χ = 1 and (b) $χ = − 1$. It can be inferred from Fig. 3 that the droplet interface appears both vertically and horizontally symmetric. In Fig. 3(a), when χ = 1, for large values of λ, which corresponds to a viscous droplet surrounded by a much less viscous fluid, the droplet is slightly squashed vertically. For small values of λ, corresponding to a viscous droplet surrounded by a much more viscous fluid, the droplet appears vertically stretched. In Fig. 3(b), when $χ = − 1$, for large values of λ, which corresponds to a viscous droplet surrounded by a much less viscous fluid, the droplet is slightly squashed vertically while for small values of λ, which corresponds to a less viscous droplet surrounded by a much more viscous fluid, the droplet appears vertically squashed in the middle. We notice that the interfacial shape in the vicinity of the north and south poles resembles the shape of a jet in the south pole of a rising bubble; the jet grows as the bubble rises and eventually collapses.
FIG. 3.

Steady-state droplet interface profiles for $ε = 1 16$ and $Ca = 1 4$ using the first 25 terms in Eq. (73) with (a) χ = 1 and (b) $χ = − 1$ for λ equal to 0.001 (red), 0.2 (blue), 0.5 (green), 0.9 (orange), 1.5 (yellow), 2.5 (brown), 5 (pink), 10 (gray), and 500 (purple).

FIG. 3.

Steady-state droplet interface profiles for $ε = 1 16$ and $Ca = 1 4$ using the first 25 terms in Eq. (73) with (a) χ = 1 and (b) $χ = − 1$ for λ equal to 0.001 (red), 0.2 (blue), 0.5 (green), 0.9 (orange), 1.5 (yellow), 2.5 (brown), 5 (pink), 10 (gray), and 500 (purple).

Close modal

Some steady-state streamlines are illustrated in Fig. 4 for $ε = 1 16 , Ca = 1 4$, and $λ = 0.5$ with (a) χ = 1 and (b) $χ = − 1$.

FIG. 4.

Steady-state streamlines for $ε = 1 16$, $Ca = 1 4$, and $λ = 0.5$ using the first 25 terms in Eq. (73) with (a) χ = 1 and (b) $χ = − 1$.

FIG. 4.

Steady-state streamlines for $ε = 1 16$, $Ca = 1 4$, and $λ = 0.5$ using the first 25 terms in Eq. (73) with (a) χ = 1 and (b) $χ = − 1$.

Close modal
We note that
$∫ − 1 1 f d η = O ( Ca 2 ),$
which is consistent with the condition $∫ − 1 1 f d η = 0$, so the volume of the droplet remains the same for all parameters.

In  Appendix E, we demonstrate that the steady-state droplet shape function f converges at various points but diverges when $η = ± 1$. In other words, the solution is valid almost everywhere except for the poles.

By differentiating Eq. (66) with respect to t and substituting in Eq. (65), we obtain
$d c j d t = − a j + 1 for j ≥ 1.$
(74)
Differentiating Eq. (67) with respect to t and subbing in Eqs. (65) and (74), we obtain
$1 Ca ( j ( j + 3 ) a j + 1 − ( j − 2 ) ( j + 1 ) a j − 1 ) = q j − 1 a d a j − 1 d t − q j + 1 a d a j + 1 d t − a j − 2 ( q j − 1 b − ξ j − 1 2 j − 1 ) − B 1 0 d 2 ζ d t 2 δ 2 , j + a j ( q j − 1 b − q j − 1 c + q j + 1 b − ξ j − 1 2 j − 1 + ξ j + 2 2 j + 3 ) + a j + 2 ( q j + 1 c − q j + 1 b − ξ j + 2 2 j + 3 ) for j ≥ 1.$
(75)
We shall now consider the solution in the small-capillary-number limit by defining
$a 1 = d K d t and a j = Ca a j ( 1 ) + O ( Ca 2 ) for j ≥ 2,$
where K is an unknown function introduced so that a1 is defined. Using this expansion, the O(1) terms in Eq. (75) are
$a j + 1 ( 1 ) = ( j − 2 ) ( j + 1 ) j ( j + 3 ) a j − 1 ( 1 ) − δ 2 , j 10 d M d t + [ δ 1 , j 4 ( q 2 b + 3 5 ξ ) − δ 3 , j 18 ( q 2 b − 2 5 ξ ) ] d K d t$
(76)
for $j ≥ 1$ where
$M = B 1 0 d ζ d t − q 1 a d K d t .$
(77)
Using Eq. (76), we obtain the recursion relation
$a 2 ( 1 ) = 1 4 ( q 2 b + 3 5 ξ ) d K d t , a 3 ( 1 ) = − 1 10 d M d t , a 4 ( 1 ) = ξ 18 d K d t , and a j + 1 ( 1 ) = ( j − 2 ) ( j + 1 ) j ( j + 3 ) a j − 1 ( 1 ) for j ≥ 4 ,$
which can be solved to yield
$a 2 j − 1 ( 1 ) = − 1 2 ( j − 1 ) ( 2 j + 1 ) d M d t for j ≥ 2 , a 2 j ( 1 ) = ξ 2 ( j + 1 ) ( 2 j − 1 ) d K d t for j ≥ 2.$
We can integrate Eq. (74) to yield
$c 1 = Ca 4 ( χ − K ) ( q 2 b + 3 5 ξ ) + O ( Ca 2 ) , c 2 j = M Ca 2 j ( 2 j + 3 ) + O ( Ca 2 ) for j ≥ 1 , c 2 j − 1 = ξ Ca ( χ − K ) 2 ( j + 1 ) ( 2 j − 1 ) + O ( Ca 2 ) for j ≥ 2.$
We can integrate Eq. (65) to yield
$b 1 = Ca 4 ( K − χ ) ( q 2 b + 3 5 ξ ) + O ( Ca 2 ) , b 2 = χ − K − Ca M 10 + O ( Ca 2 ) , b 3 = Ca 4 ( χ − K ) ( q 2 b + 17 ξ 45 ) + O ( Ca 2 ) , b 2 j = Ca M ( 4 j + 1 ) 2 ( j − 1 ) j ( 2 j + 1 ) ( 2 j + 3 ) + O ( Ca 2 ) for j ≥ 2 , b 2 j − 1 = Ca ξ ( χ − K ) ( 4 j − 1 ) 2 j ( j + 1 ) ( 2 j − 3 ) ( 2 j − 1 ) + O ( Ca 2 ) for j ≥ 3.$
Using Eq. (68), we have
$C d = 4 π B 1 0 + 4 ε π B 1 0 d K d t + Ca ( K − χ ) 6 ε π ( q 2 b + 3 5 ξ ) ( 105 − 155 B 1 0 + 58 ( B 1 0 ) 2 ) .$
We notice that the condition for the volume of the droplet to remain $4 π / 3$ given by Eq. (69) can be written as
$∑ i = 1 ∞ b 2 i + 1 = − b 1 .$
Substituting this into Eq. (66) when j = 1 gives
$b 1 + c 1 = 0 ,$
which is already satisfied. Using Eq. (59), we find that to leading order
$Λ 1 = Ca 4 ( K − χ ) ( q 2 b + 3 5 ξ ) + O ( Ca 2 ) , Λ 2 = 3 χ − 3 K − 3 Ca M 10 + O ( Ca 2 ) , Λ 3 = 3 Ca 2 ( χ − K ) ( q 2 b + 17 ξ 45 ) + O ( Ca 2 ) , Λ 2 j = Ca M ( 4 j + 1 ) 2 ( j − 1 ) ( 2 j + 3 ) + O ( Ca 2 ) for j ≥ 2 , Λ 2 j − 1 = Ca ξ ( χ − K ) ( 4 j − 1 ) 2 ( j + 1 ) ( 2 j − 3 ) + O ( Ca 2 ) for j ≥ 3.$
Using Eq. (70), we find that the vertical coordinate of the center of mass of the droplet is
$z c = 3 ε ( χ − K ) + O ( Ca 2 ) .$
Using Eq. (42), we have
$f = Ca 4 ( K − χ ) ( q 2 b + 3 5 ξ ) + ( χ − K − Ca M 10 ) 3 Q 2 Q 1 + 3 Ca 2 ( χ − K ) ( q 2 b + 17 ξ 45 ) Q 3 Q 1 + ∑ j = 2 ∞ Ca M ( 4 j + 1 ) 2 ( j − 1 ) ( 2 j + 3 ) Q 2 j Q 1 + ∑ j = 3 ∞ Ca ξ ( χ − K ) ( 4 j − 1 ) 2 ( j + 1 ) ( 2 j − 3 ) Q 2 j − 1 Q 1 + O ( Ca 2 ) ,$
(78)
where the shape of the droplet is given by $r = 1 + ε f$. Using Eq. (C12), we find that at $η = ± 1$, we have
$f | η = ± 1 = − 5 Ca 4 ( K − χ ) ( q 2 b + ξ 3 ) ± 3 ( χ − K − Ca M 10 ) ± ∑ j = 2 ∞ Ca M ( 4 j + 1 ) 2 ( j − 1 ) ( 2 j + 3 ) + ∑ j = 3 ∞ Ca ξ ( χ − K ) ( 4 j − 1 ) 2 ( j + 1 ) ( 2 j − 3 )$
plus terms of $O ( Ca 2 )$. Hence, again our solution is not valid at the poles. Notice our solution still involves two unknowns: K and the constant χ. To proceed, we shall consider two cases.
If we assume that the shape of the droplet only includes the odd modes in η, then $K = χ$, so that $a 1 = O ( Ca )$. This means that $M = B 1 0 d ζ d t$ and the solution is now defined. We notice
$a 2 j ( 1 ) = 0 for j ≥ 1 , a 3 ( 1 ) = − B 1 0 10 d 2 ζ d t 2 , a 2 j − 1 ( 1 ) = − B 1 0 2 ( j − 1 ) ( 2 j + 1 ) d 2 ζ d t 2 for j ≥ 2.$
Then, $c 2 j − 1 = O ( Ca 2 )$ for $j ≥ 1$ along with
$c 2 j = B 1 0 Ca 2 j ( 2 j + 3 ) d ζ d t + O ( Ca 2 ) for j ≥ 1.$
Then, $b 2 j − 1 = O ( Ca 2 )$ for $j ≥ 1$ along with
$b 2 = − Ca B 1 0 10 d ζ d t + O ( Ca 2 ) , b 2 j = Ca B 1 0 ( 4 j + 1 ) 2 ( j − 1 ) j ( 2 j + 1 ) ( 2 j + 3 ) d ζ d t + O ( Ca 2 )$
for $j ≥ 2$. This leads to $C d = 4 π B 1 0 + O ( Ca 2 )$. We find that
$Λ 2 j − 1 = O ( Ca 2 ) for j ≥ 1 , Λ 2 = − 3 Ca B 1 0 10 d ζ d t + O ( Ca 2 ) , Λ 2 j = Ca B 1 0 ( 4 j + 1 ) 2 ( j − 1 ) ( 2 j + 3 ) d ζ d t + O ( Ca 2 ) for j ≥ 2.$
Using Eq. (70), we find that the vertical coordinate of the center of mass of the droplet is $z c = O ( Ca 2 )$, so the center of mass of the droplet does not change in our moving coordinate system. Using Eq. (42), we then have
$f = ( − 3 Q 2 5 + ∑ j = 2 ∞ ( 4 j + 1 ) Q 2 j ( j − 1 ) ( 2 j + 3 ) ) Ca B 1 0 2 Q 1 d ζ d t$
(79)
plus terms of $O ( Ca 2 )$. Figure 5 illustrates the droplet shape using $h = ε R c sin ( π t ) / ( A c Fr 2 )$, $λ = 3 2 , ε = 0.03$, and $Ca = 1 4$ at various times. Figure 5 appears to shows that the droplet interface resembles three-dimensional axisymmetric heart-shaped solutions oscillating vertically in time, when only odd modes are present. Furthermore, in Fig. 6, we display the streaklines for the same parameters at times $1 2$, 1, and $3 2$. The streaklines are contours of $ψ ̂ 0 + ε ψ ̂ 1$ for r < 1 and contours of $ψ 0 + ε ψ 1$ for r > 1 using Eqs. (34), (35), (44), and (45). Figure 6 shows the flow pattern within the droplet and reveal as expected that the center of the recirculation zone within the droplet appears to move to the vertical position where the droplet is widest.
FIG. 5.

Unsteady droplet interfaces for case (i) using the first 25 terms in Eq. (79) for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 3 2 , ε = 0.03$, and $Ca = 1 4$ at times $1 2$ (red), $3 4$ (blue), 0.9 (green), 1 (orange), 1.1 (yellow), $5 4$ (brown), and $3 2$ (pink).

FIG. 5.

Unsteady droplet interfaces for case (i) using the first 25 terms in Eq. (79) for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 3 2 , ε = 0.03$, and $Ca = 1 4$ at times $1 2$ (red), $3 4$ (blue), 0.9 (green), 1 (orange), 1.1 (yellow), $5 4$ (brown), and $3 2$ (pink).

Close modal
FIG. 6.

Streaklines for case (i) using the first 25 terms in the series for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 3 2 , ε = 0.03$ and $Ca = 1 4$ at times (a) $1 2$, (b) 1, and (c) $3 2$.

FIG. 6.

Streaklines for case (i) using the first 25 terms in the series for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 3 2 , ε = 0.03$ and $Ca = 1 4$ at times (a) $1 2$, (b) 1, and (c) $3 2$.

Close modal
To remove the odd modes from the infinite sum in Eq. (78), we require M = 0; however, the single $Q 2 / Q 1$ odd term will remain. Hence, $B 1 0 d ζ d t = q 1 a d K d t$. Notice that $q 1 a = 3 B 1 0$, so $K = ζ / 3$. The solution is now defined up to an integration constant χ. Using Eq. (76), we obtain
$a 2 j − 1 ( 1 ) = 0 for j ≥ 2 , a 2 ( 1 ) = ( q 2 b + 3 5 ξ ) 1 12 d ζ d t , a 2 j ( 1 ) = ξ 6 ( j + 1 ) ( 2 j − 1 ) d ζ d t for j ≥ 2.$
Then, $c 2 j = O ( Ca 2 )$ for $j ≥ 1$ along with
$c 1 = Ca 12 ( 3 χ − ζ ) ( q 2 b + 3 5 ξ ) + O ( Ca 2 ) , c 2 j − 1 = ξ Ca ( 3 χ − ζ ) 6 ( j + 1 ) ( 2 j − 1 ) + O ( Ca 2 ) for j ≥ 2.$
We can integrate Eq. (65) to yield
$b 1 = Ca 12 ( ζ − 3 χ ) ( q 2 b + 3 5 ξ ) + O ( Ca 2 ) , b 2 = χ − ζ 3 + O ( Ca 2 ) , b 3 = Ca 12 ( 3 χ − ζ ) ( q 2 b + 17 ξ 45 ) + O ( Ca 2 ) , b 2 j = O ( Ca 2 ) for j ≥ 2 , b 2 j − 1 = Ca ξ ( 3 χ − ζ ) ( 4 j − 1 ) 6 j ( j + 1 ) ( 2 j − 3 ) ( 2 j − 1 ) + O ( Ca 2 ) for j ≥ 3.$
Using Eq. (68), we have
$C d = 4 π B 1 0 + 4 ε π B 1 0 3 d ζ d t + Ca ( ζ − 3 χ ) 18 ε π ( q 2 b + 3 5 ξ ) ( 105 − 155 B 1 0 + 58 ( B 1 0 ) 2 ) .$
Using Eq. (59), we find that to leading order
$Λ 1 = Ca 12 ( ζ − 3 χ ) ( q 2 b + 3 5 ξ ) + O ( Ca 2 ) Λ 2 = 3 χ − ζ + O ( Ca 2 ) , Λ 3 = Ca 2 ( 3 χ − ζ ) ( q 2 b + 17 ξ 45 ) + O ( Ca 2 ) Λ 2 j = O ( Ca 2 ) for j ≥ 2 Λ 2 j − 1 = Ca ξ ( 3 χ − ζ ) ( 4 j − 1 ) 6 ( j + 1 ) ( 2 j − 3 ) + O ( Ca 2 ) for j ≥ 3.$
Using Eq. (70), we find that the vertical coordinate of the center of mass of the droplet is $z c = ε ( 3 χ − ζ ) + O ( Ca 2 )$. Using Eq. (42) and removing terms $O ( Ca 2 )$, we get
$f = ( ζ − 3 χ ) Ca 12 ( q 2 b + 3 ξ 5 ) − ( ζ − 3 χ ) Q 2 Q 1 − ( ζ − 3 χ ) Ca 6 Q 1 ( 3 q 2 b + 17 ξ 15 ) Q 3 − ( ζ − 3 χ ) Ca 6 Q 1 ∑ j = 3 ∞ ξ ( 4 j − 1 ) Q 2 j − 1 ( j + 1 ) ( 2 j − 3 ) .$
(80)
It is evident from Eq. (80) that the value of χ, without loss of generality, can be chosen to be zero or absorbed into ζ. Figure 7 depicts the droplet shape using $h = ε R c sin ( π t ) / ( A c Fr 2 )$, $λ = 3 2$, χ = 0, $ε = 0.03$, and $Ca = 1 4$ at various times.
FIG. 7.

Unsteady-state droplet interfaces for case (ii) using the first 25 terms in Eq. (80) for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 3 2$, χ = 0, $ε = 0.03$ and $Ca = 1 4$ at times 0 (red), $1 4$ (blue), 0.4 (green), $1 2$ (orange), 0.6 (yellow), $3 4$ (brown), and 1 (pink).

FIG. 7.

Unsteady-state droplet interfaces for case (ii) using the first 25 terms in Eq. (80) for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 3 2$, χ = 0, $ε = 0.03$ and $Ca = 1 4$ at times 0 (red), $1 4$ (blue), 0.4 (green), $1 2$ (orange), 0.6 (yellow), $3 4$ (brown), and 1 (pink).

Close modal

Figure 7 shows that the droplet undergoes axisymmetric stretching and squeezing in time, when only even modes are present. Finally, Fig. 8 plots the streaklines $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 1 2 , ε = 0.01$, and $Ca = 1 4$ at times $1 2$, 1, and $3 2$. Figure 8 shows the flow pattern remains almost symmetrical throughout the oscillations.

FIG. 8.

Streaklines for case (ii) using the first 25 terms in the series for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 1 2 , ε = 0.01$, and $Ca = 1 4$ at times (a) $1 2$, (b) 1, and (c) $3 2$.

FIG. 8.

Streaklines for case (ii) using the first 25 terms in the series for $h = ε R c sin ( π t ) / ( A c Fr 2 )$, with $λ = 1 2 , ε = 0.01$, and $Ca = 1 4$ at times (a) $1 2$, (b) 1, and (c) $3 2$.

Close modal

In this study, we have considered the effect of a vertically oscillating flow field on a viscous slightly deformable axisymmetric droplet falling through a fluid using the Stokes flow equations. An expansion in the amplitude of the deformation of the droplet allowed the equations to be expanded. The zeroth-order solution yielded the well-known solutions by Rybczynski17 and Hadamard.18 The first-order equations we obtained for a droplet in a vertically oscillating flow field in axisymmetric spherical polar coordinates by expressing the solution to Stokes flow as a series involving the modified Gegenbauer polynomials. Using the interfacial conditions along with the far-field conditions and the velocities remain finite at the origin, the problem was reduced to an infinite system of linear ordinary differential equations. In the small-capillary limit, the system was analytically solved, which led to a droplet with singularities at its poles. Additionally, the drag in the vertical direction and center of mass of the droplet was obtained.

Three-dimensional axisymmetric vertically oscillating heart-shaped solutions were obtained when only odd modes are present. Heart-shaped solutions of droplets have been found in the literature, see Sostarecz and Belmonte,32 Norouzi and Davoodi.37 In the case when only even modes are present, the droplet now exhibits axisymmetric stretching and squeezing. As experiments have shown that droplets do deform from a sphere, our results do appear to be consistent with the early stages of the deformation of a spherical droplet. Experiments and numerical simulations by Pozrikidis22 and Machu et al.23 found that droplets deformations are present and increase in time. Such results may explain why the solution exhibited singularities at the poles of the droplet. One notes that singularities are common in fluid flows; for example, cusps were found by Joseph et al.,40 when partially submerged cylinders are rotated, and by Jeong and Moffatt,41 when fully submerged cylinders are rotated. Hence, the breakdown of the solution at the poles should have been expected.

This study was not able to obtain the conditions for when the droplet resonates as only the first-order terms in the expansion were obtained. Perhaps, if the expansion had been taken further, such a condition would have emerged. It would be interesting to see whether experiments could provide additional insights into the early stages of how a spherical droplet deforms in a vertically oscillating flow field. In order for the experiments to agree with this study, they would need to ensure that a pair of immiscible fluids are chosen that satisfy Re =  $O ( ε 2 )$, Fr =  $O ( ε − 1 )$, and $A c / R c = O ( ε − 1 )$ along with a small capillary number. One expects that such an experiment would require a relatively small droplet falling through a very viscous fluid. If the droplet radius was too large, one would expect much larger variations in the droplet shape than considered in this study. Furthermore, if the liquid was not sufficiently viscous, the use of Stokes flow would not be appropriate. If the density ratio is too large, i.e., when κ is large, then inertial effects would come into play, which have been neglected in this study.

The authors would like to thank Professor Damien Foster for fruitful discussions and the referees for improving the manuscript, EC Horizon 2020 Framework Programme (H2020): RISE-ATM2BT-824022.

The authors have no conflicts to disclose.

Iolo T. Williams: Formal analysis (equal); Methodology (equal). Serafim Kalliadasis: Writing – review & editing (equal). Sotos Constantinou Generalis: Writing – review & editing (equal). Philip M.J. Trevelyan: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

We write the velocities in terms of derivatives of the streamfunctions ψ and $ψ ̂$ using
$u = ∇ × A , u ̂ = ∇ × A ̂ ,$
where
$A = ψ ( r , θ ) r sin ( θ ) e ϕ , A ̂ = ψ ̂ ( r , θ ) r sin ( θ ) e ϕ$
and $e ϕ$ is a unit vector perpendicular to the plane of motion (azimuthal direction). These satisfy $∇ · u = 0$ and $∇ · u ̂ = 0$. The velocity components are given by
$u r = ψ θ r 2 sin ( θ ) , u θ = − ψ r r sin ( θ ) ,$
(A1)
$u ̂ r = ψ ̂ θ r 2 sin ( θ ) , u ̂ θ = − ψ ̂ r r sin ( θ ) .$
(A2)
We find that
$∇ × u = − ∇ 2 A , ∇ 2 u = ∇ × ( ∇ 2 A ) , ∇ × u ̂ = − ∇ 2 A ̂ , ∇ 2 u ̂ = ∇ × ( ∇ 2 A ̂ ) .$
Hence, the Stokes flow equations become
$∇ × ( ∇ 2 A ) = ∇ p and λ ∇ × ( ∇ 2 A ̂ ) = ∇ p ̂,$
and we obtain
$∇ 4 A = 0 and ∇ 4 A ̂ = 0 .$
Now,
$∇ 2 A = 1 r ( ∂ ∂ r ( ψ r sin ( θ ) ) + ∂ ∂ θ ( ψ θ r 2 sin ( θ ) ) ) e ϕ , ∇ 2 A ̂ = 1 r ( ∂ ∂ r ( ψ ̂ r sin ( θ ) ) + ∂ ∂ θ ( ψ ̂ θ r 2 sin ( θ ) ) ) e ϕ,$
and so, by introducing the operator E2 defined by
$E 2 ψ = sin ( θ ) ( ∂ ∂ r ( ψ r sin ( θ ) ) + ∂ ∂ θ ( ψ θ r 2 sin ( θ ) ) ) ,$
we can write
$∇ 2 A = E 2 ψ r sin ( θ ) e ϕ , ∇ 4 A = E 4 ψ r sin ( θ ) e ϕ , ∇ 2 A ̂ = E 2 ψ ̂ r sin ( θ ) e ϕ , ∇ 4 A ̂ = E 4 ψ ̂ r sin ( θ ) e ϕ .$
Hence, we have
$E 4 ψ = 0 and E 4 ψ ̂ = 0.$
(A3)
Finally, we introduce $η = cos ( θ )$ to get
$E 2 ψ = ∂ 2 ψ ∂ r 2 + 1 − η 2 r 2 ∂ 2 ψ ∂ η 2 .$
To solve $E 4 ψ = 0$, we let $w = E 2 ψ$ so that $E 2 w = 0$. To solve the later partial differential equation, we use the method of separation of variables and let $w = R j ( r ) Q j ( η )$, and the equation is rewritten as
$R j ″ Q j + 1 − η 2 r 2 R j Q j ″ = 0 ⇒ r 2 R j ″ R j = ( η 2 − 1 ) Q j ″ Q j = − I .$
The separation constant is given by $I = − j ( j + 1 )$, where j is a non-negative integer, and we find that
$R j = A ̃ r j + 1 + B ̃ r − j , ( η 2 − 1 ) Q j ″ − j ( j + 1 ) Q j = 0.$
(B1)
By letting $L j = Q j ′$, we see that
$( 1 − η 2 ) L j ″ − 2 η L j ′ + j ( j + 1 ) L j = 0,$
which is the Legendre differential equation. The Ljs are orthogonal functions with a weight function of 1 [see Leal39 Eq. (7-120)], since
$∫ − 1 1 L i L j d η = 2 δ i , j 2 j + 1,$
where Kronecker delta $δ i , j$ is zero unless i = j in which case it becomes 1. The Ljs are the Legendre polynomials. The first few Legendre polynomials are as follows:
$L 0 = 1 , L 1 = η , L 2 = 3 2 η 2 − 1 2 , L 3 = 5 2 η 3 − 3 2 η , …$
We note that we are assuming the flow to be axisymmetric so along with no $ϕ$ dependence, we have $u θ = 0$ on θ = 0 and $θ = π$, i.e., at $η = ± 1$. This implies $ψ r = 0$ at $η = ± 1$; hence, ψ is constant at $η = ± 1$. From the far-field condition, the constant vanishes, hence,
$ψ = 0 at η = ± 1.$
We therefore define
$Q j = ∫ − 1 η L j ( η ) d η$
so that Qj = 0 at $η = − 1$. Here, Qj are a modified set of orthogonal Gegenbauer polynomials. The first few such polynomials are
$Q 0 = η + 1 , Q 1 = 1 2 ( η 2 − 1 ) , Q 2 Q 1 = η , Q 3 Q 1 = 1 4 ( 5 η 2 − 1 ) , …$
Hence, the general solution for w is given by
$w = ∑ j = 1 ∞ ( A ̃ j r j + 1 + B ̃ j r − j ) Q j .$
It should be noted that we start at j = 1, since Q0 does not satisfy the axisymmetric condition. To solve now solve $w = E 2 ψ$, we let $ψ = ψ I + ψ c$ and
$E 2 ψ I = w = ( A ̃ j r j + 1 + B ̃ j r − j ) Q j .$
(B2)
At the same time, we notice that
$E 2 ( Ξ r σ Q j ) = Ξ r σ − 2 ( σ ( σ − 1 ) − j ( j + 1 ) ) Q j ,$
using Eq. (B1). Now, we compare the powers of r. If we let $Ξ j = A j$ with $j + 1 = σ − 2 ⇒ σ = j + 3$, then $A ̃ j = A j ( 4 j + 6 )$. If we let $Ξ j = B j$ with $− j = σ − 2 ⇒ σ = 2 − j$, then $B ̃ j = B j ( 2 − 4 j )$. Hence,
$ψ I = ∑ j = 1 ∞ ( A j r j + 3 + B j r 2 − j ) Q j ,$
where
$E 2 ψ I = ∑ j = 1 ∞ ( A j ( 4 j + 6 ) r j + 1 + B j ( 2 − 4 j ) r − j ) Q j .$
Now, $E 2 ψ c = 0$, so that
$ψ c = ∑ j = 1 ∞ ( C j r j + 1 + D j r − j ) Q j .$
Hence, we have found that
$ψ = ∑ j = 1 ∞ ( A j r j + 3 + C j r j + 1 + B j r 2 − j + D j r − j ) Q j ,$
(B3)
$ψ ̂ = ∑ j = 1 ∞ ( A ̂ j r j + 3 + C ̂ j r j + 1 + B ̂ j r 2 − j + D ̂ j r − j ) Q j .$
(B4)
Now, $∇ 2 u = ∇ p$, so $∇ × ( ∇ × u ) = − ∇ p$. Hence,
$∇ × u = − 1 r sin ( θ ) ( E 2 ψ ) e ϕ,$
which implies that
$∇ p = ∇ × ( 1 r sin ( θ ) ( E 2 ψ ) e ϕ ) .$
By using the curl in spherical polar coordinates, we obtain
$∇ p = 1 r 2 sin ( θ ) ∂ ( E 2 ψ ) ∂ θ e r − 1 r sin ( θ ) ∂ ( E 2 ψ ) ∂ r e θ .$
By using the gradient in spherical polar coordinates, we obtain
$∂ p ∂ r = 1 r 2 sin ( θ ) ∂ ( E 2 ψ ) ∂ θ and 1 r ∂ p ∂ θ = − 1 r sin ( θ ) ∂ ( E 2 ψ ) ∂ r .$
Using $η = cos ( θ )$, we obtain
$∂ p ∂ r = − 1 r 2 ∂ ( E 2 ψ ) ∂ η and ∂ p ∂ η = 1 1 − η 2 ∂ ( E 2 ψ ) ∂ r .$
Now, $w = E 2 ψ I$, where
$E 2 ψ I = ∑ j = 1 ∞ ( A j ( 4 j + 6 ) r j + 1 + B j ( 2 − 4 j ) r − j ) Q j .$
Solving for p yields
$p = p 0 + ∑ j = 1 ∞ ( − A j 4 j + 6 j r j + B j 2 − 4 j j + 1 r − j − 1 ) Q j ′,$
(B5)
and similarly,
$p ̂ = p ̂ 0 + λ ∑ j = 1 ∞ ( − A ̂ j 4 j + 6 j r j + B ̂ j 2 − 4 j j + 1 r − j − 1 ) Q j ′ .$
(B6)
The Qjs are a modified set of Gegenbauer polynomials satisfying the ordinary differential equation
$2 Q 1 Q j ″ = j ( j + 1 ) Q j for j ∈ N$
(C1)
along with Qj = 0 at $η = ± 1$. The Qjs are orthogonal functions with a weight function of $Q 1 − 1$ [Leal39 Eq. (7-123)] since
$∫ − 1 1 Q i Q j Q 1 d η = − 4 δ i , j j ( j + 1 ) ( 2 j + 1 ) .$
Similarly, $Q j ′$ s and the $Q j ″$ s are linearly independent when applying the weight functions 1 and Q1, respectively. We find that
$∫ − 1 1 Q i ′ Q j ′ d η = 2 δ i , j 2 j + 1 , ∫ − 1 1 Q 1 Q i ″ Q j ″ d η = − j ( j + 1 ) δ i , j 2 j + 1 .$
We notice that the polynomials satisfy the recursion relations
$Q j + 1 = j Q 1 ′ Q j + 2 Q 1 Q j ′ j + 2 for j ≥ 0,$
(C2)
$Q j − 1 = ( j + 1 ) Q 1 ′ Q j − 2 Q 1 Q j ′ j − 1 for j ≥ 2.$
(C3)
For $j ≥ 1$, we have
$( Q 1 2 ( Q j Q 1 ) ′ ) ′ ≡ Q j ″ Q 1 − Q j = 1 2 ( j + 2 ) ( j − 1 ) Q j ,$
(C4)
$Q 1 Q j ′ = ( j + 1 ) ( j + 2 ) Q j + 1 − j ( j − 1 ) Q j − 1 2 ( 2 j + 1 ) ,$
(C5)
$Q 1 ′ Q j = ( j + 2 ) Q j + 1 + ( j − 1 ) Q j − 1 2 j + 1 .$
(C6)
Furthermore, for $j ≥ 3$, we have
$Q 1 Q 1 ′ ( Q j Q 1 ) ′ = ∑ m = 1 j − 2 1 + ( − 1 ) j − m 2 j ( j + 1 ) m ( m + 1 ) ( 2 m + 1 ) Q m + ( j − 1 ) Q j .$
(C7)
We notice that
$∫ − 1 1 Q j d η = − 2 3 δ 1 , j , ∫ − 1 1 η Q j ′ d η = 2 3 δ 1 , j ,$
(C8)
$∫ − 1 1 Q 2 j Q 1 d η = 0 , ∫ − 1 1 Q 2 j − 1 Q 1 d η = 2 j ( 2 j − 1 ) ,$
(C9)
$∫ − 1 1 Q 1 ′ Q 2 j − 1 Q 1 d η = 0 , ∫ − 1 1 Q 1 ′ Q 2 j Q 1 d η = 2 j ( 2 j + 1 ) .$
(C10)
Furthermore, we notice that
$Q k Q 1 = 1 on η = 1 ,$
(C11)
$Q k Q 1 = ( − 1 ) k + 1 on η = − 1.$
(C12)
We can write each Qj as a sum in powers of η, i.e.,
$Q j = ∑ m = 0 j + 1 k m η m,$
where
$k m + 2 = k m m ( m − 1 ) − j ( j + 1 ) ( m + 1 ) ( m + 2 ) for 0 ≤ m ≤ j − 1.$
If j is even, then $k 0 = 0$ and $k 1 = ( − 1 ) j / 2 ( j − 1 ) ! 2 j − 1 ( j 2 ) ! ( j 2 − 1 ) !$, while if j is odd, then $k 1 = 0$ and $k 0 = ( − 1 ) ( j + 1 ) / 2 ( j − 1 ) ! 2 j ( j + 1 2 ) ! ( j − 1 2 ) !$.
The equation for the surface of the droplet is given by Eq. (27), and so, we have
$r θ = ( 1 + ε f ) e θ − ε f η sin ( θ ) e r and r ϕ = ( 1 + ε f ) sin ( θ ) e ϕ .$
The vector cross product of these vectors is
$r θ × r ϕ = ( 1 + ε f ) 2 sin ( θ ) e r + ε f η ( 1 + ε f ) sin 2 ( θ ) e θ .$
Hence,
$| r θ × r ϕ | = ( 1 + ε f ) sin ( θ ) ( 1 + ε f ) 2 + ε 2 f η 2 sin 2 ( θ ) ,$
which reduces to $| r θ × r ϕ | ≈ ( 1 + 2 ε f ) sin ( θ )$ in the small ε-limit. We then have
$C d = ∫ 0 2 π ∫ 0 π τ ¯ n · e z ( 1 + 2 ε f ) sin ( θ ) d θ d ϕ = 2 π ∫ − 1 1 τ ¯ n · e z ( 1 + 2 ε f ) d η = 2 π ∫ − 1 1 1 + 2 ε f N η ( 2 e r r − p ) d η + 4 π ∫ − 1 1 1 + 2 ε f N ( ε r η f η − 1 ) 1 − η 2 e r θ d η − 2 π ∫ − 1 1 1 + 2 ε f N ( ε r f η ( 1 − η 2 ) ( 2 e θ θ − p ) ) d η = 2 π ∫ − 1 1 1 + 2 ε f N η [ 4 r 3 ψ η − 2 ψ r η r 2 − p ] d η + 4 π ∫ − 1 1 1 + 2 ε f r N [ ε r η f η − 1 ] [ ψ r r − ψ r r 2 − Q 1 r 2 ψ η η ] d η + 8 π ε ∫ − 1 1 1 + 2 ε f r 3 N f η Q 1 ( ψ r η − η ψ r 2 Q 1 − ψ η r − p 2 r 2 ) d η ,$
where r is evaluated on the surface of the droplet, i.e., Eq. (19), which is substituted in and linearizing we obtain
$C d = 2 π ∫ − 1 1 η ( 4 ψ η 0 − 2 ψ r η 0 − p 0 ) − 2 ψ r 0 + ψ r r 0 + 2 Q 1 ψ η η 0 d η + 2 ε π ∫ − 1 1 η ( 4 ψ η 1 − 2 ψ r η 1 − p 1 ) − 2 ψ r 1 + ψ r r 1 + 2 Q 1 ψ η η 1 d η + 2 ε π ∫ − 1 1 f η ( 4 ψ η r 0 − 4 ψ η 0 − 2 ψ r r η 0 − p r 0 − 2 p 0 ) d η − 2 ε π ∫ − 1 1 f η η [ ψ r r 0 + 2 Q 1 ψ η η 0 ] d η − 4 ε π ∫ − 1 1 f η Q 1 ( p 0 + 2 ψ η 0 − 2 ψ r η 0 ) d η + 2 ε π ∫ − 1 1 f [ − ψ r r 0 + ψ rrr 0 + 2 Q 1 ( ψ η η r 0 − ψ η η 0 ) ] d η .$
(D1)
Using integration by parts on the last line in Eq. (D1) and simplifying, we obtain
$C d = 2 π ∫ − 1 1 η ( 4 ψ η 0 − 2 ψ r η 0 − p 0 ) − 2 ψ r 0 + ψ r r 0 + 2 Q 1 ψ η η 0 d η + 2 ε π ∫ − 1 1 η ( 4 ψ η 1 − 2 ψ r η 1 − p 1 ) − 2 ψ r 1 + ψ r r 1 + 2 Q 1 ψ η η 1 d η + 2 ε π ∫ − 1 1 f [ 2 Q 1 p η 0 − η p r 0 + ( 4 η 2 − 2 ) ψ η η 0 ] d η + 2 ε π ∫ − 1 1 f [ ψ rrr 0 − η ψ r r η 0 − 2 Q 1 ψ η η r 0 ] d η .$
Using the zeroth-order solutions in Eqs. (34), (35), and (38), we obtain
$C d = π ∫ − 1 1 η 2 ( 9 − 6 B 1 0 ) − η p 0 + 3 B 1 0 − 3 d η + 2 ε π ∫ − 1 1 η ( 4 ψ η 1 − 2 ψ r η 1 − p 1 ) − 2 ψ r 1 + ψ r r 1 + 2 Q 1 ψ η η 1 d η,$
which further reduces to
$C d = 4 π B 1 0 + 2 ε π ∫ − 1 1 η ( 4 ψ η 1 − 2 ψ r η 1 − p 1 ) d η − 2 ε π ∫ − 1 1 2 ψ r 1 − ψ r r 1 − 2 Q 1 ψ η η 1 d η .$
Substituting in the first-order solutions (44), (45), and (53) and simplifying yields
$C d = 4 π B 1 0 + 4 ε π ∫ − 1 1 η ∑ j = 1 ∞ Q j ′ ( A ̂ j 1 + C ̂ j 1 ) [ j − 2 + 3 j + 1 ] d η + 4 ε π ∫ − 1 1 η ∑ j = 1 ∞ Q j ′ D j 1 [ 4 − 3 j + 1 ] d η + 2 ε π ∫ − 1 1 ∑ j = 1 ∞ Q j ( ( A ̂ j 1 + C ̂ j 1 ) ( 2 j 2 − 2 ) + D j 1 ( 2 + 4 j ) ) d η .$
Finally, using the equations in Eq. (C8), we obtain
$C d = 4 π B 1 0 + 4 3 ε π ( A ̂ 1 1 + C ̂ 1 1 − D 1 1 ) = 4 π B 1 0 + 4 3 ε π ( 3 B 1 0 a 1 + α β c 1 ) + 2 3 ε π ( 15 − 23 B 1 0 + 10 ( B 1 0 ) 2 ) b 1 .$
(D2)
Evaluating Eq. (73) at η = 0 yields
$f | η = 0 = − χ Ca 4 ( ξ + 5 2 q 2 b ) + O ( Ca 2 )$
using the value of k0 in  Appendix C. This shows that f converges at η = 0. The solution appears to converge for other values of η as well, for example, at $η = 1 2$,
$f | η = 1 2 ≈ 3 χ 2 − 5 χ Ca q 2 b 32 − 0.162 847 χ Ca ξ + O ( Ca 2 )$
and at $η = cos ( π 12 ) = 6 + 2 4 ≈ 0.965 926$,
$f | η = 6 + 2 4 ≈ 3 χ 6 + 2 4 + 5 χ Ca q 2 b 8 ( 5 + 3 3 ) + 2.991 395 χ Ca ξ + O ( Ca 2 ) .$
However, using Eq. (C12), we find that
$f | η = ± 1 = 3 χ ± 5 χ Ca 12 ( 3 q 2 b + ξ ) ± ∑ j = 3 ∞ χ Ca ξ ( 4 j − 1 ) 2 ( j + 1 ) ( 2 j − 3 ) + O ( Ca 2 ) .$
Notice each term in the series is of $O ( j − 1 )$ for large j, so the series does not converge, which means that f does not exist at $η = ± 1$.
1.
D.
Chatterjee
,
P.
Jain
, and
K.
Sarkar
, “
Ultrasound-mediated destruction of contrast microbubbles used for medical imaging and drug delivery
,”
Phys. Fluids
17
,
100603
(
2005
).
2.
M.
Ward
,
J.
Wu
, and
J.
Chiu
, “
Ultrasound-induced cell lysis and sonoporation enhanced by contrast agents
,”
J. Acoust. Soc. Am.
105
,
2951
2957
(
1999
).
3.
J.
Krehbiel
,
L.
Schideman
,
D.
King
, and
J.
Freund
, “
Algal cell disruption using microbubbles to localize ultrasonic energy
,”
Bioresour. Technol.
173
,
448
451
(
2014
).
4.
G. G.
Stokes
, “
On the effect of the internal friction of fluids on the motion of pendulums
,”
Trans. Cambridge Philos. Soc.
9
,
8
106
(
1851
).
5.
C. W.
Oseen
, “
Über den gültigkeitsbereich der stokesschen widerstandsformel
,”
Ark. Mat. Astron. Fys.
9
,
1
16
(
1913
).
6.
I.
Proudman
and
J.
Pearson
, “
Expansions at small Reynolds number for the flow past a sphere and a circular cylinder
,”
J. Fluid Mech.
2
,
237
262
(
1957
).
7.
L. E.
Payne
and
W. H.
Pell
, “
The Stokes flow problem for a class of axially symmetric bodies
,”
J. Fluid Mech.
7
,
529
549
(
1960
).
8.
R.
Cox
, “
The steady motion of a particle of arbitrary shape at small Reynolds numbers
,”
J. Fluid Mech.
23
,
625
643
(
1965
).
9.
J. R.
Ockendon
, “
The unsteady motion of a small sphere in a viscous liquid
,”
J. Fluid Mech.
34
,
229
239
(
1968
).
10.
W.
Chester
,
D.
Breach
, and
I.
Proudman
, “
On the flow past a sphere at low Reynolds number
,”
J. Fluid Mech.
37
,
751
760
(
1969
).
11.
H.
Pruppacher
,
B.
Clair
, and
A.
Hamielec
, “
Some relations between drag and flow pattern of viscous flow past a sphere and a cylinder at low and intermediate Reynolds numbers
,”
J. Fluid Mech.
44
,
781
790
(
1970
).
12.
L.
Landau
and
E.
Lifshitz
,
Course of Theoretical Physics
, Fluid Mechanics Vol.
6
(
Pergamon
,
New York
,
1987
).
13.
R.
Mei
and
R.
, “
Flow past a sphere with an oscillation in the free-stream velocity and unsteady drag at finite Reynolds number
,”
J. Fluid Mech.
237
,
323
341
(
1992
).
14.
C.
Chang
,
B.-H.
Liou
, and
R.-L.
Chern
, “
An analytical and numerical study of axisymmetric flow around spheroids
,”
J. Fluid Mech.
234
,
219
246
(
1992
).
15.
H.
Taşeli
and
M.
Demiralp
, “
A new approach to the classical stokes flow problem: Part I Methodology and first-order analytical results
,”
J. Comput. Appl. Math.
78
,
231
232
(
1997
).
16.
S.
Otto
, “
On stability of the flow around an oscillating sphere
,”
J. Fluid Mech.
239
,
47
68
(
1992
).
17.
W.
Rybczynski
, “
Über die fortschreitende bewegung einer flüssigen kugel in einem zähen medium
,”
A
,
40
46
(
1911
).
18.
J. S.
, “
Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquid visqueux
,”
152
,
1735
1738
(
1911
).
19.
T. D.
Taylor
and
A.
Acrivos
, “
On the deformation and drag of a falling viscous drop at low Reynolds number
,”
J. Fluid Mech.
18
,
466
476
(
1964
).
20.
S.
Lin
and
A.
Gautesen
, “
Creeping flow around a deforming sphere
,”
J. Fluid Mech.
56
,
61
71
(
1972
).
21.
D.
Oliver
and
J.
Chung
, “
Steady flows inside and around a fluid sphere at low Reynolds numbers
,”
J. Fluid Mech.
154
,
215
230
(
1985
).
22.
C.
Pozrikidis
, “
The instability of a moving viscous drop
,”
J. Fluid Mech.
210
,
1
21
(
1990
).
23.
G.
Machu
,
W.
Meile
,
L. C.
Nitsche
, and
U.
Schaflinger
, “
Coalescence, torus formation, and breakup of sedimenting drops: Experiments and computer simulations
,”
J. Fluid Mech.
447
,
299
336
(
2001
).
24.
D. K.
Srivastava
,
R. R.
, and
S.
, “
Steady stokes flow around deformed sphere. Class of oblate axisymmetric bodies
,”
Int. J. Appl. Math. Mech.
8
,
17
53
(
2012
).
25.
J. D.
Krehbiel
and
J. B.
Freund
, “
Stokes flow inside a sphere in an inviscid extensional flow
,”
Z. Angew. Math. Phys.
68
,
1
13
(
2017
).
26.
S.
Sahu
and
A. S.
Khair
, “
Dynamics of a viscous drop under an oscillatory uniaxial extensional stokes flow
,”
Int. J. Multiphase Flow
146
,
103844
(
2022
).
27.
H.
Godé
,
S.
Charton
,
E.
Climent
, and
D.
Legendre
, “
Basset-Boussinesq history force acting on a drop in an oscillatory flow
,”
Phys. Rev. Fluids
8
,
073605
(
2023
).
28.
F. M.
Leslie
and
R. I.
Tanner
, “
The slow flow of a visco-elastic liquid past a sphere
,”
Q. J. Mech. Appl. Math.
14
,
36
48
(
1961
).
29.
B.
Caswell
and
W. H.
Schwarz
, “
The creeping motion of a non-Newtonian fluid past a sphere
,”
J. Fluid Mech.
13
,
417
426
(
1962
).
30.
A.
Beris
,
J.
Tsamopoulos
,
R.
Armstrong
, and
R.
Brown
, “
Creeping motion of a sphere through a Bingham plastic
,”
J. Fluid Mech.
158
,
219
244
(
1985
).
31.
H.
Ramkissoon
, “
Stokes flow past a non-Newtonian fluid spheroid
,”
Z. Angew. Math. Mech.
78
,
61
66
(
1998
).
32.
M. C.
Sostarecz
and
A.
Belmonte
, “
Motion and shape of a viscoelastic drop falling through a viscous fluid
,”
J. Fluid Mech.
497
,
235
252
(
2003
).
33.
R. B.
Bird
,
R. C.
Armstrong
, and
O.
Hassager
,
Dynamics of Polymeric Liquids
, Fluid Mechanics Vol.
1
, 2nd ed. (
Wiley-Interscience
,
1987
).
34.
S.
Mukherjee
and
K.
Sarkar
, “
Viscoelastic drop falling through a viscous medium
,”
Phys. Fluids
23
,
013101
(
2011
).
35.
B. R.
Jaiswal
and
B.
Gupta
, “
Wall effects on Reiner-Rivlin liquid spheroid
,”
Appl. Comput. Mech.
8
,
157
176
(
2014
).
36.
B.
Vamerzani
,
M.
Norouzi
, and
B.
, “
Theoretical and experimental study on the motion and shape of viscoelastic falling drops through Newtonian media
,”
Rheol. Acta
55
,
935
955
(
2016
).
37.
M.
Norouzi
and
M.
Davoodi
, “
Analytical study on motion and shape of creeping Boger drops falling through viscoelastic media
,”
J. Braz. Soc. Mech. Sci. Eng.
40
,
125
(
2018
).
38.
B. R.
Jaiswal
, “
Steady stokes flow of a non-Newtonian Reiner-Rivlin fluid streaming over an approximate liquid spheroid
,”
Appl. Comput. Mech.
14
,
145
162
(
2020
).
39.
L.
Leal
,
(
Cambridge University Press
,
Cambridge
,
2007
).
40.
D. D.
Joseph
,
J.
Nelson
,
M.
Renardy
, and
Y.
Renardy
, “
Two-dimensional cusped interfaces
,”
J. Fluid Mech.
223
,
383
409
(
1991
).
41.
J.-T.
Jeong
and
H. K.
Moffatt
, “
Free-surface cusps associated with flow at low Reynolds number
,”
J. Fluid Mech.
241
,
1
22
(
1992
).
Published open access through an agreement withJISC Collections