A falling fluid droplet in an oscillating flow field

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 axi-symmetric stretching and squeezing.


I. INTRODUCTION
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, Stokes 4 examined steady flow past a solid sphere of radius R c and moving at uniform speed U in the absence of inertia and he obtained the stream function of the flow field w ¼ ðUR c =4rÞð3r 2 À R 2 c Þ sin 2 ðhÞ, where r is the radial coordinate measured from the center of the sphere and h 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 ¼ 6plR c U, where l is the dynamic viscosity of the fluid.Stokes 4 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.
Oseen 5 obtained the first-order correction to flow past a solid sphere for low Reynolds numbers, with the dimensionless streamfunction given by and, further, that the drag coefficient on the sphere in the vertical direction was Proudman and Pearson 6 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 ¼ 6p Re ð1 þ 3Re 8 þ 9 40 Re 2 lnðReÞ þ OðRe 2 ÞÞ.Payne and Pell 7 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.Cox 8 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.Ockendon 9 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 ¼ 6p Re ð1 þ 3Re 8 þ 9 40 Re 2 ½lnðReÞ þ c þ 27 80 Re 3 lnðReÞ þ OðRe 2 ÞÞ, where c ¼ c þ 5 3 lnð2Þ À 323 360 and c 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 Lifshitz 12 gave the solution for an oscillating spherical drop in an infinite medium and found the smallest possible frequency of oscillations of the drop was ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi , where a is the surface tension coefficient, q is the density of the fluid, and R c 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 Adrian 13 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 Demiralp 15 examined axisymmetric Stokes flow past an arbitrary axisymmetrical solid body by writing the solution as an infinite series involving Gegenbauer polynomials.Otto 16 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€ ortler vortices.
The flow of an fluid sphere through a another fluid in the absence of inertia was analyzed by Rybczynski 17 and Hadamard 18 who independently found that the dimensionless streamfunction inside ŵ and outside w of the sphere were given by where B ¼ ð2l þ 3lÞ=ð2l þ 2lÞ, where l and l 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 pðq À qÞgR 3 c where R c is the radius of the sphere, g is the magnitude of the gravitational acceleration, and q and q denote the density of the fluids inside and outside the sphere, respectively.Furthermore, the speed of the sphere was given by U ¼ gR 2 c ðq À qÞ=ð3BlÞ.Taylor and Acrivos 19 theoretically investigated the axisymmetric motion of a slightly deformable fluid drop falling through a fluid in the small-butfinite 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 Gautesen 20 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 Chung 21 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.Pozrikidis 22 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 Freund 25 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 Khair 26 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 e 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 Tanner 28 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 Schwarz 29 looked at low-Reynolds-number flow of a nondeformable 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.Ramkissoon 31 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 Belmonte 32 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 Sarkar 34 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 Gupta 35 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 Davoodi 37 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.Jaiswala 38 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 smallcapillary limit.Finally, a summary of our findings and conclusions are offered in Sec.X.

II. MODEL EQUATIONS
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 where hats denote fluid 2 (the droplet), U is the fluid velocity, T is time, q 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.
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, A c 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 where A c dH=dT is the time-dependent forcing velocity and 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., 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 =DT ¼ 0, where D/DT is the convective derivative, or The stress balance equation on the interface is where s ¼ 2lẽ is the deviatoric stress tensor, with ẽ ¼ 1 2 ð$U þ ð$UÞ T Þ the rate-of-strain tensor, c is the surface tension, n is the unit outward pointing normal vector to the interface, and l ¼ q is the dynamic viscosity.We shall assume that the surface tension is constant, and so, the tangential and normal stress balances are ððs À ŝ ÞnÞ Á t i ¼ 0 on F ¼ 0; (6) 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: where lower case letters without tildes are used to denote the dimensionless parameters, R c is the average radius of the droplet, and U c 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ð2pU 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 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 We have assumed that j is not too large and k 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 where . Now, we use Fðx; y; z; tÞ ¼ 0 as the equation for the interface.The continuity of velocity condition equation ( 4) becomes and the kinematic condition equation ( 5) becomes Finally, the tangential stress condition equation ( 6) is and the normal stress condition equation ( 7) as where where Ca is the capillary number and Fr is the Froude number.

III. AXISYMMETRIC SPHERICAL POLAR COORDINATES
To solve the Stokes flow equations in Eq. ( 10), we shall use axisymmetric spherical polar coordinates and write the velocities as where e r and e h are the unit vectors in the r and h directions, respectively.We recall that r is the distance measured from the origin and h is the angle measured anticlockwise from the positive z axis.For the functions ûr ðr; hÞ and ûh ðr; hÞ we require, ûr and ûh are bounded at r ¼ 0: (16)   Using e z ¼ cosðhÞe r À sinðhÞe h , the far-field condition (11) becomes u r !cosðhÞ and u h !ÀsinðhÞ as r !1; and the continuity of velocity condition (12) becomes We now suppose that where e 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 and the kinematic condition ( 13) is written as Now, the unit normal and tangential vectors to the interface, F ¼ 0, are where Finally, the tangential stress condition ( 14) is written as where the components of the rate-of-strain tensor are as follows: @u r @h : The normal stress condition ( 15) is given by By introducing g ¼ cosðhÞ, we can express the velocity components in terms of the streamfunctions w and ŵ for the flow fields outside and inside of the droplet, respectively, as The Stokes flow equations (10) are satisfied when the streamfunctions satisfy equation (A3), namely, where which is explained in the book by Leal. 39The boundedness condition at r ¼ 0 for the velocity components, Eq. ( 16), means that ŵ The far-field condition (17) becomes By integrating these expressions and using We now let f ðh; tÞ ¼ f ðg; tÞ, which makes the equation for the interface in vector form.The continuity of velocity condition (18) becomes Expanding each term using a Taylor series about e ¼ 0 and linearizing in e, we obtain The kinematic condition (20) now becomes Expanding about e ¼ 0 and linearizing in e yield The tangential stress condition ( 21) is also linearized in e, and we obtain We can express the components of the rate-of-strain tensor in terms of the stream function as The tangential stress condition is written as Upon expanding about e ¼ 0 and linearizing in e, we get on r ¼ 1.We now turn to the normal stress condition (22), which is also linearized in e to give which upon substituting in the components of the rate-of-strain tensor yields 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ðe À1 Þ; Fr ¼ Oðe À1 Þ, and Re ¼ Oðe 2 Þ, with f and h both of O(1).Then, expanding about e ¼ 0 and linearizing, we obtain To solve the above system of equations, we introduce the expansions wðr; g; tÞ ¼ w 0 ðr; gÞ þ ew

IV. ZEROTH-ORDER SOLUTION
We begin by seeking the zeroth-order solutions, which are independent of time and the droplet is a sphere of radius 1.The zerothorder Stokes flow equations are whose general solutions are given by Eqs.(B3) and (B4), namely, where the Q j s 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 B0 j ¼ D0 j ¼ 0 for j ! 1.The far-field condition (26) then yields This means that A 0 j ¼ 0 for j !1; C 0 1 ¼ À1, and C 0 j ¼ 0 for j ! 2. The kinematic condition (30) yields As the derivatives of the Gegenbauer polynomials are linearly independent, we can equate their coefficients to yield Next, the continuity of velocity conditions ( 28) and (29) gives The first condition gives We next turn to the tangential stress condition (31), which is written as and B 0 j ¼ 0 for j ! 2. At this stage, the leading-order stream functions are given by Since g ¼ cosðhÞ and Q 1 ¼ 1 2 ðg 2 À 1Þ, we can write 2 , and so, w 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= ffiffi ffi 2 p , with the minimum value ŵ0 ¼ Àb=16.This means that a recirculation zone exists inside the droplet.We notice that if k ! 1, then B 0 1 !3=2 and ŵ0 !0. Figure 2 illustrates the streamlines obtained using Eqs.(34)   and (35).The general solutions for p 0 and p0 are given by Eqs.(B5) and (B6), namely, and as presented in Appendix B and shown by in the book Leal. 39ubstituting in the values obtained for the coefficients A 0 j ; B 0 j ; Â0 j , and B0 j into Eqs.( 36) and (37) yields The normal stress condition (32) yields Substituting in the solutions for p0 ; p 0 , w 0 , and ŵ0 into Eq.( 39) and equating coefficients of the constants and Finally, using the definitions of the dimensionless parameters, the second equation above becomes Substituting in the ratios for j and k yields 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

V. FIRST-ORDER EQUATIONS
We next turn our attention to the first-order solutions in which the droplet interface is allowed to deform in g and time.The firstorder Stokes flow equations are whose general solutions are The boundedness condition for the velocities at r ¼ 0 (25) yields ŵ1 r 2 is bounded at r ¼ 0; which means that B1 j ¼ D1 j ¼ 0 for j ! 1.The far-field condition (26) yields Turning then to the continuity of velocity conditions ( 28) and (29), Substituting in w 0 and ŵ0 reduces these to where a ¼ 12B 0 1 À 15.As the Q j s are linearly independent functions, the first condition implies that where for j ! 1.Here, we are assuming that a 6 ¼ 0, i.e., B 0 1 6 ¼ 5 4 , i.e., k 6 ¼ 1, i.e., l 6 ¼ l and the two fluids have different dynamic viscosities.We notice that at this stage, we can write the first-order streamfunctions as Turning now to the kinematic condition (30), we have Substituting in w 0 reduces this to from which where for j ! 1.We now make use of the identity (C5) in Appendix C to obtain By defining X 0 ¼ 0 and since the Q j s are linearly independent, we obtain We now turn to the tangential stress condition (31), and substituting in w 0 and ŵ0 leads to which becomes where we have defined Next, we use the shape of the droplet, Eq. ( 42), along with the identity (C7) in Appendix C to get where we have defined as for j ! 1.By changing the order of the summations and using the fact that the Q j s are linearly independent functions, the tangential stress condition reduces to where for j ! 1.By setting m ¼ 2i þ j, the above expression is simplified to We next look at the normal stress condition (32) and substituting in w 0 and ŵ0 , from Eqs. ( 34) and ( 35) and p 0 and p0 from Eq. ( 38), and using Eq. ( 40) lead to where n ¼ 30 À 21B 0 1 and f ¼ AcFr 2 eRc dh dt .Using the general solutions for the pressure, we have Substituting the pressures into our normal stress condition yields where we have defined as We now use the shape of the droplet in Eq. (42) to get 2 Ca To simplify this expression, we make use of the identities (C4)-(C7) in Appendix C to give 2 Ca As before, we change the order of the summations and simplify the above equation to 1 Ca where we have defined K 0 ¼ H 0 ¼ 0. As the Q j s are linearly independent, the normal stress condition reduces to where we have substituted (50) in the tangential stress condition, and d i;j is the Kronecker delta, equal to 1 when i ¼ j and zero otherwise.

VI. SIMPLIFYING THE LINEAR SYSTEM
Equations ( 43), ( 46), (49), and (54) are expressions for K j , X j , C j , and H j in terms of Â1 j ; Ĉ1 j , and D 1 j .We notice that we can invert Eqs.(43), (46), and (49) to yield To express the system of Eqs. ( 47), ( 50), (51), and (55) in a simpler form, we define the new variables where a 0 ¼ b 0 ¼ d 0 ¼ 0. Furthermore, using the new variables, we can write where (63) for j ! 1 with q a 0 ¼ q b 0 ¼ q c 0 ¼ 0. The kinematic condition equation (47) becomes Substituting Eq. ( 51) into the tangential stress condition equation ( 50) with these new variables yields Next, the normal stress condition equation (55 for j ! 1.We can eliminate d j from this equation to obtain Hence, our system of equations consists now of Eqs. ( 65)-(67) which we need to solve to obtain a j , b j, and c j .Once b j is known, we can obtain K j from Eq. ( 59), and then, using Eq. ( 42), we can obtain f, and then, from r ¼ 1 þ ef , we can obtain the shape of the droplet.Using Eqs. ( 56)-( 58), we can obtain Â1 j ; Ĉ1 j , and D 1 j .The streamfunctions and the pressures can then be obtained from Eqs. ( 44), (45), and (53).

VII. DRAG, VOLUME, AND CENTER OF MASS OF DROPLET
The drag on the droplet only acts in the z direction and is given by where the drag coefficient is given by In Appendix D, we find that the drag coefficient is given by Eq. ( 68), namely, Notice, the relation for the drag coefficient reduces to the one for the unforced case that of a free droplet given by Leal. 39The volume of the droplet is given by Linearizing in e gives Using Eqs. ( 42) and (C9), we obtain Conservation of mass dictates that for the droplet to have a volume of 4p=3, corresponding to the volume of the undisturbed droplet, a sphere of radius 1, we require that 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 z c is given by so that the droplet volume remains as 4p=3.Using Eq. (C10), we obtain for the vertical coordinate of the center of the droplet.

VIII. SMALL CAPILLARY NUMBER-STEADY-STATE SOLUTION
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 This means that a 2j ¼ 0 and a 2jÀ1 ¼ a 1 for j ! 1.However, if the a 2jÀ1 s are all equal, then this means that the X 2jÀ1 s tend to infinity as j tends to infinity, which means that at least one of the Âj s, Ĉj s, or D j s would diverge.Hence, in order to obtain a convergent solution, we require that a 1 ¼ 0. Equation (66) reduces to Substituting this into Eq.( 67) yields Using Eq. ( 69) for a steady state, we have For convenience, we introduce the constant I 0 defined as Then, using the odd values in Eq. ( 71) and simplifying we have for j ! 1. Next, using the even values in Eq. ( 71) and simplifying we have for j ! 1.To solve these in the small-capillary-number limit, we let b 2 ¼ v, where v 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 Substituting these into the odd equations at leading order gives 2jð2j À 3Þb Solving this yields b for j !3: Using now Eq. ( 66), we have Hence, Using Eq. ( 68) for the drag coefficient, we have We notice that when v > 0, the OðCaÞ term is reducing the drag coefficient as long as % 57:05.Using Eq. ( 59), we obtain Using Eq. ( 70), we find that the vertical coordinate of the center of mass of the droplet is given by z c ¼ 3ev þ OðCa 2 Þ.Furthermore, from Eq. ( 42), we obtain Using x ¼ r cosðhÞ and z ¼ r sinðhÞ with g ¼ cosðhÞ and r ¼ 1 þ ef , 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 e ¼ 1 16 and Ca ¼ 1 4 for various values of k with (a) v ¼ 1 and (b) v ¼ À1.It can be inferred from Fig. 3 that the droplet interface appears both vertically and horizontally symmetric.In Fig. 3(a), when v ¼ 1, for large values of k, which corresponds to a viscous droplet surrounded by a much less viscous fluid, the droplet is slightly squashed vertically.For small values of k, corresponding to a viscous droplet surrounded by a much more viscous fluid, the droplet appears vertically stretched.In Fig. 3(b), when v ¼ À1, for large values of k, which corresponds to a viscous droplet surrounded by a much less viscous fluid, the droplet is slightly squashed vertically while for small values of k, 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.Some steady-state streamlines are illustrated in Fig. 4 for e ¼ We note that ð 1 which is consistent with the condition Ð 1 À1 fdg ¼ 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 g ¼ 61.In other words, the solution is valid almost everywhere except for the poles.
Differentiating Eq. ( 67) with respect to t and subbing in Eqs. ( 65) and (74), we obtain We shall now consider the solution in the small-capillary-number limit by defining where K is an unknown function introduced so that a 1 is defined.Using this expansion, the O(1) terms in Eq. ( 75) are for j ! 1 where Using Eq. (76), we obtain the recursion relation and a which can be solved to yield We can integrate Eq. (74) to yield We can integrate Eq. (65) to yield Using Eq. ( 68), we have We notice that the condition for the volume of the droplet to remain 4p=3 given by Eq. ( 69) can be written as Substituting this into Eq.(66) when j ¼ 1 gives which is already satisfied.Using Eq. ( 59), we find that to leading order Using Eq. (70), we find that the vertical coordinate of the center of mass of the droplet is Using Eq. ( 42), we have where the shape of the droplet is given by r ¼ 1 þ ef .Using Eq. (C12), we find that at g ¼ 61, we have 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 v.To proceed, we shall consider two cases.

A. Case (i): Odd modes
If we assume that the shape of the droplet only includes the odd modes in g, then K ¼ v, so that a 1 ¼ OðCaÞ.This means that M ¼ B 0 1 df dt and the solution is now defined.We notice 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 plus terms of OðCa 2 Þ. Figure 5 illustrates the droplet shape using h ¼ eR c sinðptÞ=ðA c Fr 2 Þ, k ¼ 3 2 ; e ¼ 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 þ e ŵ1 for r < 1 and contours of w 0 þ ew 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.
Using Eq. ( 68), we have Using Eq. ( 59), we find that to leading order Using Eq. ( 70), we find that the vertical coordinate of the center of mass of the droplet is z c ¼ eð3v À fÞ þ OðCa 2 Þ.Using Eq. ( 42) and removing terms OðCa 2 Þ, we get It is evident from Eq. ( 80) that the value of v, without loss of generality, can be chosen to be zero or absorbed into f. Figure 7 depicts the droplet shape using h ¼ eR c sinðptÞ=ðA c Fr 2 Þ, k ¼ 3 2 , v ¼ 0, e ¼ 0:03, and Ca ¼ 1 4 at various times.
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 ¼ eR c sinðptÞ=ðA c Fr 2 Þ, with k ¼ 1 2 ; e ¼ 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.

X. DISCUSSION AND CONCLUSIONS
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 Rybczynski 17 and Hadamard. 18The firstorder 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 heartshaped 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. 37In 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 Pozrikidis 22 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ðe 2 Þ, Fr ¼ Oðe À1 Þ, and A c =R c ¼ Oðe À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 j is large, then inertial effects would come into play, which have been neglected in this study.
w h r 2 sinðhÞ e / ; and so, by introducing the operator E 2 defined by we can write r sinðhÞ e / : Hence, we have Finally, we introduce g ¼ cosðhÞ to get To solve E 4 w ¼ 0, we let w ¼ E 2 w 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 ðgÞ, and the equation is rewritten as The separation constant is given by I ¼ Àjðj þ 1Þ, where j is a nonnegative integer, and we find that which is the Legendre differential equation.The L j s are orthogonal functions with a weight function of 1 [see Leal 39 Eq. (7-120)], since where Kronecker delta d i;j is zero unless i ¼ j in which case it becomes 1.The L j s are the Legendre polynomials.The first few Legendre polynomials are as follows: We note that we are assuming the flow to be axisymmetric so along with no / dependence, we have u h ¼ 0 on h ¼ 0 and h ¼ p, i.e., at g ¼ 61.This implies w r ¼ 0 at g ¼ 61; hence, w is constant at g ¼ 61.From the far-field condition, the constant vanishes, hence, w ¼ 0 at g ¼ 61: We therefore define Here, Q j are a modified set of orthogonal Gegenbauer polynomials.The first few such polynomials are Hence, the general solution for w is given by It should be noted that we start at j ¼ 1, since Q 0 does not satisfy the axisymmetric condition.To solve now solve w ¼ E 2 w, we let w ¼ w I þ w c and At the same time, we notice that using Eq.(B1).Now, we compare the powers of r.If we let N where Now, E 2 w c ¼ 0, so that Hence, we have found that which implies that By using the curl in spherical polar coordinates, we obtain By using the gradient in spherical polar coordinates, we obtain Using g ¼ cosðhÞ, we obtain

APPENDIX C: PROPERTIES OF THE GEGENBAUER POLYNOMIALS
The Q j s are a modified set of Gegenbauer polynomials satisfying the ordinary differential equation along with Q j ¼ 0 at g ¼ 61.The Q j s are orthogonal functions with a weight function of Q À1 1 [Leal 39 Eq. (7-123)] since ð 1 Similarly, Q 0 j s and the Q 00 j s are linearly independent when applying the weight functions 1 and Q 1 , respectively.We find that ð 1 We notice that the polynomials satisfy the recursion relations For j ! 1, we have Furthermore, for j ! 3, we have We notice that Furthermore, we notice that We can write each Q j as a sum in powers of g, i.e., where for 0 m j À 1: If j is even, then k 0 ¼ 0 and k 1 ¼ ðÀ1Þ j=2 ðjÀ1Þ!

APPENDIX D: DRAG ON THE DROPLET
The equation for the surface of the droplet is given by Eq. ( 27), and so, we have The vector cross product of these vectors is q ; which reduces to jr h Â r / j % ð1 þ 2ef Þ sinðhÞ in the small e-limit.
We then have where r is evaluated on the surface of the droplet, i.e., Eq. ( 19), which is substituted in and linearizing we obtain f g 4w 0 gr À 4w 0 g À 2w 0 rrg À p 0 r À 2p 0 dg Using integration by parts on the last line in Eq. (D1) and simplifying, we obtain f ½w 0 rrr À gw 0 rrg À 2Q 1 w 0 ggr dg: Using the zeroth-order solutions in Eqs.(34), (35), and (38), we obtain which further reduces to Substituting in the first-order solutions (44), (45), and (53) and simplifying yields Finally, using the equations in Eq. (C8), we obtain

APPENDIX E: EVALUATING THE STEADY-STATE DROPLET SHAPE FUNCTION f
Evaluating Eq. (73) at g ¼ 0 yields using the value of k 0 in Appendix C.This shows that f converges at g ¼ 0. The solution appears to converge for other values of g as well, for example, at g ¼ However, using Eq.(C12), we find that 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 g ¼ 61.

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