PHYSICAL REVIEW E 90, 023007 (2014)

Stability analysis of a thin liquid film on an axially oscillating cylindrical surface in the high-frequency limit Selin Duruk and Alexander Oron* Department of Mechanical Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel (Received 17 February 2014; revised manuscript received 3 July 2014; published 13 August 2014) We consider an axisymmetric liquid film on a horizontal cylindrical surface subjected to axial harmonic oscillation in the high-frequency limit. We derive and analyze the nonlinear evolution equation describing the nonlinear dynamics of this physical system in terms of the averaged film thickness. The method used for the derivation of the evolution equation is based on long-wave theory and the separation of the relevant fields into fast and slow components. We carry out the linear stability analysis for a film of a constant thickness which shows that axial forcing of the cylinder may result in either stabilization or destabilization of the axisymmetric flow with respect to the unforced one, depending on the choice of the parameter set. The analysis is extended to the weakly nonlinear stage and it reveals that the system bifurcates subcritically from the equilibrium. DOI: 10.1103/PhysRevE.90.023007

PACS number(s): 47.20.Dr, 47.20.Ma

I. INTRODUCTION

External forcing is very well known to introduce a dramatic change in the behavior of mechanical systems. The best example of such striking influence is the inverted pendulum [1,2], which is also known as Kapitza’s inverted pendulum. Recently, several intriguing experimental observations concerning the behavior of liquid systems, both films and drops, were reported in the literature [3–5]. Investigation of the influence of a planar substrate oscillation on the dynamics of a horizontal liquid layer overlying it dates more than four decades back with the linear stability analysis of Yih [6], who carried out long-wave expansions to determine domains of instability. Later, Or [7] revisited Yih’s problem and demonstrated, using Floquet theory, that the instability sets in first for a finite-wavelength regime in a wide range of parameters. Halpern and Frenkel [8] studied the nonlinear stability of a horizontal bilayer system subjected to the Rayleigh-Taylor instability when the substrate on the top of the heavier fluid layer performs in-plane oscillations. The dynamics of the interface separating the two layers was investigated in terms of a weakly nonlinear evolution equation of the Kuramoto-Sivashinsky type derived there for both low-frequency and high-frequency limits. In certain parameter domains, the Rayleigh-Taylor instability was shown to saturate due to planar oscillations. Investigation of linear stability of thin liquid films on a cylindrical surface was carried out by Goren [9]. He derived a dispersion relation in terms of Bessel functions and determined the fastest growing modes in the general case as well as in the limits of negligible inertia and negligible fluid viscosity. He also favorably compared the results of linear stability analysis with his own experimental observations. Hammond [10] considered the nonlinear evolution of axisymmetric films on the inner wall of a horizontal cylinder. He showed that the film is prone to capillary instability and breaks up in the

*

[email protected]

1539-3755/2014/90(2)/023007(13)

long-time limit into a series of almost disconnected liquid droplets which due to their axial symmetry represent collars on a cylinder. Gauglitz and Radke [11] extended Hammond’s theory by introducing a more accurate treatment to the interfacial curvature term. Yarin, Oron, and Rosenau [12] and Johnson et al. [13] developed nonlinear theories describing the dynamics of thin films on the outer and inner surfaces of a horizontal cylinder using the integral method based on spatial averaging in the radial direction. Weidner, Schwartz, and Eres [14] considered a three-dimensional evolution of a thin liquid film on a cylindrical surface. They identified topologically different drop structures for the case of a larger cylinder diameter where the drops tend to emerge on the underside of the cylindrical surface and for the case of a smaller cylinder diameter where the drops assumed a more symmetric configuration with a significant portion of liquid located on the upper side of the cylinder. They also characterized various stages in time evolution of the surface-tension and gravitational potential energies as well as of that for viscous dissipation. Axially uniform films on a horizontal cylindrical surface were studied by Reisfeld and Bankoff [15] who discussed possible steady states and determined the time evolution of the film in the gravity-dominated regime in terms of elliptic functions. Haimovich and Oron [16] considered the nonlinear dynamics of an axisymmetric film on a horizontal cylindrical surface under “low” -frequency harmonic axial vibration and found that the film rupture taking place on a static cylindrical surface is prevented at the nonlinear stage provided that the product of the amplitude and frequency of forcing exceeds a certain value depending only on the liquid properties and the geometry. These results were extended and generalized for the case of a horizontal cylinder subjected to differential heating with respect to the ambient gas phase [17] and to the case of a double-frequency forcing in the isothermal case [18]. Several papers were devoted to the investigation of the dynamics of films in the limit of high frequency in the cases of harmonic vibration in the directions normal, tangential, and oblique to the planar horizontal substrate and for either stabilizing or destabilizing gravity [19–23].

023007-1

©2014 American Physical Society

SELIN DURUK AND ALEXANDER ORON

PHYSICAL REVIEW E 90, 023007 (2014)

It is a purpose of this paper to develop a mathematical model describing the nonlinear dynamics of a thin axisymmetric liquid film coating a horizontal cylindrical surface undergoing harmonic axial vibration in the limit of a “high” frequency using the methods of long-wave theory [24,25] and multiple scaling similar to those employed by Shklyaev, Alabuzhev, and Khenner [20]. This investigation contains the linear and weakly nonlinear stability analyses for a uniform thickness film in the framework of the model, whereas the fully nonlinear analysis of the system will be carried out elsewhere. The paper is organized as follows. Section II is devoted to the problem statement and the formulation of the governing equations. Section III deals with the decomposition of the equations and boundary conditions into separate problems for the averaged and pulsating components and with their respective solutions. Section IV contains the formulation of the governing evolution equation and the discussion related to its various limiting cases. Section V discusses the results of linear stability analysis for a uniform-thickness base state. Section VI addresses the weakly nonlinear analysis of the problem based on the evolution equation derived in Sec. IV. Section VII presents the concluding remarks. II. PROBLEM STATEMENT AND GOVERNING EQUATIONS

We consider a thin liquid film of the average thickness h0 coating a horizontal solid circular cylinder of radius Rc . The liquid is assumed to be incompressible, isothermal, and with the following properties: density ρ, dynamic viscosity μ, kinematic viscosity ν = μ/ρ, and surface tension σ , all assumed to be constant. We assume the solid substrate to perform axial harmonic vibration with the displacement amplitude b and a single frequency ω. In what follows, we assume that the liquid film is axisymmetric, so gravity effects can be neglected in comparison with the capillary effects. This is feasible if the cylinder radius and the mean film thickness are each much smaller than the σ 1/2 capillary length of the liquid c = ( ρg ) . In this case, the 2

relevant Bond number, Bo =  2 , with the typical length size c of  = Rc + h0 , will be much smaller than one, and gravity is dominated by capillarity. This result fits the results obtained by de Bryun [26]. As pointed by Hammond [10], in the linear stage of evolution the nonaxisymmetric modes are more stable than the axisymmetric one. This, of course, does not preclude the emergence of nonaxisymmetric patterns in the nonlinear stage of the film evolution. A full and thorough account for the nonaxisymmetric effects leads to a tremendous increase in the numerical effort related to the solution of the appropriate twodimensional evolution equation. Also, the numerical results obtained by Weidner, Schwartz, and Eres [14], who studied nonaxisymmetric films on a horizontal static cylindrical surface, show that, for cylinders of a smaller diameter, patterns tend to be more symmetric with respect to the cylinder axis. The governing equations are therefore the axisymmetric version of the continuity and the Navier-Stokes equations written in cylindrical coordinates with x and r being the axial and radial spatial coordinates, respectively, where the

azimuthal component of the velocity field is zero: (rw)r + rux = 0,



ρ(wt + wwr + uwx ) = −pr + μ  ρ(ut + wur + uux ) = −px + μ

∂ ∂r



(1a)  1 ∂ (rw) + wxx , r ∂r 



1 ∂ (rur ) + uxx , r ∂r

(1b) (1c)

where w and u represent, respectively, the r− and x− components of the velocity field, t is time, and p is pressure. The boundary conditions are as follows: at the surface of the cylinder r = Rc , u = bω cos ωt, w = 0,

(2a)

and at the free surface r = h + Rc , ht − w + uhx = 0, Tij nj = −κσ ni , i,j = r,x,

(2b) (2c)

where Tij are components of the stress tensor T, ni are those of the unit outer normal vector n to the free surface, and κ is the local twice mean interfacial curvature given by, respectively,    −1/2 hxx (−1,hx ) 1 1 + h2x − n=  , κ= . 2 Rc + h 1 + h x 1 + h2x (3) Equation (2c) is now rewritten via its radial (normal to the interface for the axisymmetric case) and axial components μhx (ur + wx ) + p − 2μwr − σ κ = 0,

(4a)

μ(ur + wx ) + hx (p − 2μux − σ κ) = 0,

(4b)

respectively. The pressure and capillary terms, respectively, p and σ κ terms, are eliminated using Eqs. (4a) and (4b), similar to Refs. [14,16], that leads to the tangential projection of the stress balance boundary condition in the form   (5) 2hx wr + 1 − h2x (ur + wx ) − 2hx ux = 0. To apply the long-wave theory we assume that the average film thickness h0 is much smaller than the characteristic wavelength l of a disturbance in the x direction and introduce the film parameter

=

h0  1. l

(6)

The equations and boundary conditions (1) and (2) are normalized by scaling the variables using h2

ν r t, r ∗ = , p∗ = 02 p, 2 Rc ρν h0 ν (w,u) = ( w ∗ ,u∗ ), (x,h) = h0 (x ∗ / ,h∗ ), h0 t∗ =

(7)

where the variables with an asterisk represent dimensionless counterparts of the corresponding dimensional ones.

023007-2

STABILITY ANALYSIS OF A THIN LIQUID FILM ON . . .

PHYSICAL REVIEW E 90, 023007 (2014)

Equations (1), (2a), (2b), (4a), and (5) are thus rewritten in scaled dimensionless form where the asterisks are dropped for simplicity,    ∂ 1 ∂ 2 2 (rw)

(wt + awwr + uwx ) = −apr + a ∂r r ∂r + O( 3 ),

(ut + awur + uux ) = − px + a 2



(8a)



1 ∂ (rur ) + 2 uxx , r ∂r (8b)

a(rw)r + rux = 0,

u = B cos( −1 t),

and at r = 1 + ah, ∂ a(1 + ah)ht + ∂x



1+ah

(9a)

ru dr = 0,

(9b)

1

  aur + 2 2ahx wr − 2hx ux − ah2x ur + wx + O( 4 ) = 0, (9c)  2 awr − aur hx − p + C

2 hxx a − 1 + ah 1 + 2 h2x

 −1/2 × 1 + 2 h2x + O( 3 ) = 0,



(9d)

where  = ωh20 /ν,

To investigate the film evolution in the high-frequency limit, we separate the equations and the boundary conditions (8) and (9) into fast and slow variables similar to what was done by Shklyaev, Alabuzhev, and Khenner [20] in the planar case and solve the problem including terms of orders up to O( ) with respect to the small parameter . We now introduce the time scales to be used in the forthcoming analysis via t¯ = t, τ = −1 t,

(8c)

at r = 1, w = 0,

III. DECOMPOSITION OF THE FIELDS INTO AVERAGED AND PULSATING PARTS

B = b/ h0 , a = h0 /Rc , C = σ h0 /ρν 2 (10)

are, respectively, the dimensionless forcing frequency and amplitude, the geometrical parameter, and inverse capillary number. We assume in this paper that a = O(1), B = O(1), and C = O(1). We will comment below [following Eq. (17c)] on the assumption related to the inverse capillary number C. As noted in the beginning of the next section, the assumption   2 will embody “the high-frequency limit” considered in the present analysis.

(11a)

serving in what follows as slow and fast times, respectively, and decompose the fields of dependent variables into their averaged and pulsating components according to ¯ t¯) +

¯ t¯) + −1 p

(t¯,τ ), h = h( h(t¯,τ ), p = p( (11b) ¯ t¯) +

w = w( ¯ t¯) + w

(t¯,τ ), u = u( u(t¯,τ ). Note that the time scale τ emerges naturally in the boundary condition Eq. (9a). To allow the times t¯ and τ to be, respectively, slow and fast, namely, t¯  τ , the requirement   2 needs to be satisfied. Therefore, to consider the film evolution in the limit of high frequency, the range of  should cover values between moderately small to order one. In the decompositions given by Eq. (11b), the values with overbar represent the averaged values of the corresponding fields X that depend only on the slow time t¯, 2π 1 ¯ X ≡ X  = X dτ (12) 2π 0 (note two alternative notations for the same value which are used upon convenience), whereas the values with a tilde correspond to the scaled pulsating parts. Note that the average value of the difference X − X  for the corresponding field X with respect to the fast time τ is zero. In what follows, we retain the leading-order terms for both pulsating and averaged parts of the fields and take into account all terms excluding those of O( 2 ) order and higher. Applied to Eqs. (8) and (9) this yields with the exception of a crucially necessary term in the interfacial curvature [see the discussion below following Eq. (17c)]

p¯ r + p˜ r = 0,

(13a)

w

ur +

u

ux ) = −( p¯ x + p˜ x ) + (

uτ + a

2

a ∂ [r( u¯ r +

ur )], r ∂r

ux ) = 0, a[r( w¯ + w

)]r + r( u¯ x +

(13b) (13c)

at r = 1,

w¯ + w

= 0,

u¯ +

u = B cos τ,

(14a)

and at r = 1 + a(h¯ +

h),

1+a h+a

1+a h+a

¯ ¯



h h ∂ ∂ 2

¯ ¯ ¯

r u¯ dr + r

u dr = 0. a(1 + a h)hτ + a(1 + a h)ht¯ + a hhτ +

∂x 1 ∂x 1 We introduce new functions via the following primitives: r r

≡ u¯ d, G(r) 

u d. F¯ (r) ≡ 1

1

023007-3

(14b)

(15)

SELIN DURUK AND ALEXANDER ORON

PHYSICAL REVIEW E 90, 023007 (2014)

In terms of these functions, Eq. (14b) is rewritten as ∂ ∂

¯

¯ h¯ t¯ + a 2 

a(1 + a h) hτ + a(1 + a h) h

hτ + [F¯ (1 + a h¯ + a

h)] + [G(1 + a h¯ + a

h)] = 0. ∂x ∂x

around r = 1 + a h¯ yields up to O( 2 ) terms in Eq. (14b) Expanding F¯ and G 

r ) 

∂(

hG ∂G ∂ F¯ 2

¯ ¯ ¯

+ a(1 + a h)ht¯ + a hhτ + +a = 0. a(1 + a h)hτ + ∂x ∂x ∂x

(16)

(17a)

The boundary condition Eq. (9c) is rewritten at r = 1 + a h¯ up to O( 2 ) terms as

u¯ r +

ur + a

h

urr = 0.

(17b)

The boundary condition Eq. (9d) projected onto r = 1 + a h¯ is written with the terms up to O( 3 ) as   h a a 2

2¯ ˜ + C − hxx − a h˜ p˜ r = 0. − ( p¯ + p) −

¯ 2 1 + a h¯ (1 + a h) Before continuing, an important remark is to be made about the presence of the 2 and 3 terms in Eq. (17c). Normally, when a long-wave expansion is applied to a film flow in the planar geometry, the inverse capillary number is assumed to be large. As a result, the emerging waves are long, surface gradients are small, and the surface-tension term is retained in the leading-order approximation. The above-mentioned promotion of the capillary term to the leading order allows one to derive a mathematically well-posed evolution equation. Without this promotion, the capillary term is of a higher order and does not enter the evolution equation at the leading order. This situation renders the evolution equation to be ill-posed containing only instability terms (backward diffusion) and solutions of this equation do not saturate due to a lack of dissipation, which could have arisen from the capillary term, and display a finite-time blowup. For a film flow on a cylindrical surface the situation is different and more intricate. Here, the surface tension affects both axial and radial (still axisymmetric) deflections, represented by several terms of disparate orders in the limit of → 0 included in the expression multiplied by the inverse capillary number C in Eq. (17c). If C is of unity order and only a term of O(1) is kept, the rescaled normal-stress-balance boundary condition Eq. (17c) will contain just one contribution of surface tension that is related only to the azimuthal curvature. Neglecting other terms and carrying out the asymptotic procedures leads at the final stage to the evolution equation with no dissipation term which must be provided by the second (longitudinal) curvature component in Eq. (3). Without this component, the situation will be similar to that in the Cartesian case addressed above without the “promotion” of the capillary term. To include surface-tension effects at the lowest order possible in in Eq. (17c) with both competing components of the interfacial curvature along with the assumption of the parameter C of a unity order, we have retained in the boundary condition Eq. (17c) all terms up to O( 3 ). These include one term of O( 2 ) arising from the projection of the location of the interface onto its averaged location, and three O( 3 ) terms. One of them is shown in Eq. (17c) and represents the contribution of the longitudinal component of the interfacial curvature. Two others which are not shown in Eq. (17c) are associated with the correction due to the interfacial metrics, and the additional contribution due to the azimuthal component

(17c)

a h¯ 2

of the interfacial curvature. These two terms, namely, 2 2(1+ax h) ¯ h2 2 a 3

2 ¯

and − (1+a h) ¯ 3 , are both proportional to hx [note that h is also proportional to h¯ x , as will follow from Eq. (23) below] and fully nonlinear, thus do not contribute anything to the linear stability analysis; see Sec. V. Moreover, they do not contribute anything sizable to the forthcoming weakly nonlinear analysis either and do not alter the value of the relevant Landau constant, see Sec. VI, and to the conclusion regarding the type of bifurcation. Therefore, only the term 2 h¯ xx is retained in Eq. (17c) due to its crucial role in the analysis. An effectively similar strategy (without separation of the fields into the averaged and pulsating components) leading to the presence of the two components of the interfacial curvature was used by Hammond [10], Weidner, Schwartz, and Eres [14], Craster and Matar [27], Novbari and Oron [28], Haimovich and Oron [16–18], and others. As will be seen below, the O( 2 ) term in Eq. (17c) will vanish in the averaging procedure. A joint appearance of both components of the interfacial curvature is known to be vital to preserve a correct stability threshold [16,27,28] of the model equations as shown in Sec. V and is also supported by comparison with experiments [27–29]. The procedure to treat the boundary value problem consisting of Eqs. (13), (14a), and (17) will be as follows: these equations are first averaged with respect to the fast time τ and only terms of O( ) are retained along with that of two orders higher in the normal stress balance Eq. (17c). Then, the averaged equations are subtracted from Eqs. (13), (14a), and (17), respectively, and only O(1)-order terms are kept to obtain a set of equations for the pulsating components which are treated first. This procedure ensures that both competing terms appearing in the curvature, Eq. (3), are included.

A. Solution for the pulsating components of the fields

Following the procedure described above, we obtain a linear set of governing equations and boundary conditions that reads

023007-4

p˜ r = 0, 

uτ = −p˜ x + a(r w

)r + r

ux = 0,

(18a) a2 ∂ (r

ur ), r ∂r

(18b) (18c)

STABILITY ANALYSIS OF A THIN LIQUID FILM ON . . .

PHYSICAL REVIEW E 90, 023007 (2014)

at r = 1, w

= 0,

u = B cos τ,

(18d)

given by the following set of equations and boundary conditions: p¯ r = 0,

¯ and at r = 1 + a h,

∂G ¯

a(1 + a h) hτ + = 0, ∂x

ur = 0,

a

w

ur  + 

u

ux  = −p¯ x +

(18e) p

= 0.

1 ∂ (r

ur ) , r ∂r and its solution is sought in the form A2

uτ =

u = Re(U eiτ ),

(19)

(20)

where U = U (x,r) is an unknown function to be determined next. Hereafter Re denotes the real part of the expression contained in the corresponding expression in the brackets. Upon substitution of this ansatz into Eq. (19) we find using boundary conditions Eqs. (18d) and (18f) arising from Eqs. (14a) and (17b), respectively, an expression for

u in the τ -periodic form   I0 (αr)K1 (Q) + K0 (αr)I1 (Q) iτ

u = B Re e , (21) I0 (α)K1 (Q) + K0 (α)I1 (Q) √ √ ¯ Here In ,Kn where α ≡ A i = i/a and Q = α(1 + a h). are, respectively, the modified Bessel functions of nth order. In fact, these Bessel functions are related to the Kelvin functions of zero and first orders, and may be in principle expressed via them. However, rather than doing this, we prefer to use the present notation. Following this, w

is obtained by integration of the continuity equation (18c) in the form

a ∂ (r u¯ r ), r ∂r

a(r w) ¯ r + r u¯ x = 0,

(18f)

The problem given by Eqs. (18) is now solved. Since p

¯ we conclude is independent of r and vanishes at r = 1 + a h, that p

is identically zero, p

= 0. Equation (18b) is rewritten in terms of A2 = /a 2 as

(24a) 2

(24b) (24c)

at r = 1, w¯ = 0,

u¯ = 0,

(24d)

¯ and at r = 1 + a h, ∂ ¯ h¯ t¯ + ∂ [F¯ (1 + a h)]+a ¯ ¯

a(1 + a h) [(1 + a h) h

u] = 0, ∂x ∂x (24e) h

urr , u¯ r = −a

 p¯ = C

 a pr

− 2 h¯ xx − a

h. 1 + a h¯

(24f) (24g)

Here X  denotes the averaged value of an arbitrary function X over the fast time τ ; see Eq. (12). Also, note that p

≡ 0, see the previous subsection; therefore, the last term in Eq. (24g) vanishes. Since the averaged pressure p¯ is independent of r, see Eq. (24a), it is determined by its value at the film interface, Eq. (24g) with p˜ = 0. Equation (24b) is integrated twice in the radial direction r, ¯ first time over the domain [1 + a h,r], second time over [1,r], and the boundary conditions Eqs. (24d) and (24f) are applied. As a result, we obtain u¯ =

w

= Re(W e )h¯ x   1 − r + α [I1 (αr)K0 (α) + K1 (αr)I0 (α)] iτ Bh¯ x = e . Re 1 + a h¯ α 2 [I1 (Q)K0 (α) + K1 (Q)I0 (α)]2 (22) iτ

(B)2 2a 2 +

p¯x 2a 2



r

[F1 () + F2 ()]d

1

 (r 2 − 1) ¯ 2 ln r , − (1 + a h) 2

(25)

Substituting

u and w

into Eq. (18e) yields the expression for

h,

where the functions F1 (r) and F2 (r) are defined via the following primitives:

h = Re(H eiτ )h¯ x   B h¯ x ieiτ . (23) = Re ¯ 2 α 2 [I0 (α)K1 (Q) + K0 (α)I1 (Q)]2 (1 + a h)

r 1 ∂ |U |2 d, F1 (r) = a h¯ x f1 = r ∂x 1+a h¯ F2 (r) = a h¯ x f2 = a Re(W U ∗ ).

It is important to emphasize that in the absence of forcing, B = 0, the pulsating components of the fields,

u, w

, and

h vanish. The latter two also vanish if the averaged film thickness h¯ is constant.

(26)

The velocity component u¯ can be recast into the form

B. Solution for the components of the averaged fields

Based on the procedure described above, the problem for the averaged components of the fields is 023007-5

r (B)2 ¯ ¯ u¯ = hx [f (,h)]d 2a 1   p¯x (r 2 − 1) ¯ 2 ln r , + 2 − (1 + a h) 2a 2

(27)

SELIN DURUK AND ALEXANDER ORON

PHYSICAL REVIEW E 90, 023007 (2014)

where ¯ ≡ f1 (,h) ¯ + f2 (,h) ¯ f (,h)  ∗ I0 (α ) [M1 K1 (Q) − M3 I1 (Q)] − K0 (α ∗ ) [M2 K1 (Q) − M4 I1 (Q)] 1 Re = |N |2 QN ∗  I0 (α ∗ )K1 (−iQ) + I1 (−iQ)K0 (α ∗ ) , − αQN ∗ ¯ = I1 (Q)K0 (α) + K1 (Q)I0 (α), N (α,h) ∗

(28)

∗

M1 (r) = I1 (αr)K0 (α r), M2 (r) = I1 (αr)I0 (α r), M3 (r) = K1 (αr)K0 (α ∗ r), M4 (r) = K1 (αr)I0 (α ∗ r).

The asterisks denote complex conjugates of the correspond¯ = N (−iα,h) ¯ was used to ing values and the property N ∗ (α,h) simplify the expressions. IV. EVOLUTION EQUATION

In this section we obtain the evolution equation describing the spatiotemporal nonlinear dynamics of the averaged value of the film thickness h¯ in a slow time t¯. Various limiting cases of this evolution equation are also discussed. The expression obtained in the previous section for the velocity component u¯ and given by Eq. (27) is now substituted into Eq. (24e). This yields the key equation of the paper in the form of the pertinent evolution equation in terms of n = 1 + ah in the form  n  ∂ nnt¯ + u d + nx D(n) = 0 ∂x 1 written explicitly as nnt¯ + +





nxxx ∂ C anx Y (n) + 2 2 2 ∂x a n a

where 2 J (n) = 2





n

r 1

r

(29a)

 ¯ f (,h)d dr,

1

   1 , D(n) = − 2 Re 3 2n α N |N |2

1 2 ∂ ∂ (n )t + [nζ L1 (n)] + [L2 (n)nζ ζ ζ ] = 0, 2 ∂ζ ∂ζ

(31a)

where L1 (n) ≡ ϒ[J (n) + a 2 D(n)] +

aY (n) , n2

L2 (n) =

1 Y (n), a (31b)

and B2 . (31c) C It is important to emphasize that function L2 (n) is always positive and that the influence of external forcing is accounted for via the component of L1 (n) proportional to ϒ. Note that the dimensionless forcing amplitude B appears in Eq. (31a) only via ϒ; however, the presence of the frequency  comes in many terms in functions D(n) and J (n). It is clear that the effect of vibration on the film dynamics will be negligible for small ϒ; however, it becomes more and more pronounced with an increase in ϒ. The second component of L1 (n) containing Y (n) is related to the destabilizing effect of the azimuthal curvature of the interface. In the case of an unforced system, B = 0, Eq. (29a) reduces to Eq. (25) of Ref. [16] where only two terms are taken in the axial derivative of the interfacial curvature. Furthermore, in this case for B = 0 and in the limit of a  1, Eq. (29a) upon a proper rescaling x → a x,t → 3 Ct2 a 3 reduces to that derived by Hammond [10], namely ϒ=



B2 ∂ [nx (J (n) + a 2 D(n))] = 0, a 2 ∂x

contribution of the pulsating components of the fields

u and

h via the kinematic boundary condition. We will return to this discussion below. We now normalize Eq. (29a) using the scaling a2 (30) t¯ = 2 t, x = ζ, C and Eq. (29a) is then rewritten in the conservative form

(29b)

Y (n) = (4n4 ln n − 3n4 + 4n2 − 1)/16. In Eq. (29a) the second term is related to capillarity, whereas the third term stands for the effect of vibration. The latter is of a diffusion type; therefore, vibration introduces extra stabilization if the “diffusion coefficient” J (n) + a 2 D(n) is negative, or extra destabilization if it is positive. Note that, in the absence of forcing, B = 0, the third term of Eq. (29a) vanishes. Hence the resulting evolution equation in the unforced case will contain only the capillary term with its two competing components related to the destabilizing azimuthal and stabilizing axial curvatures. Among the terms with the factor B 2 arising from forcing, the first term is related to the convective contribution of the averaged part of the ¯ whereas the second one is associated with the flow field u,

ht + (h3 hx + h3 hxxx )x = 0.

(32)

In the case of a → 0 corresponding to the planar (Rc → ∞) configuration with forcing (B = 0), the value α has a in its denominator, so effectively the arguments of the Bessel functions involved here tend to infinity. Therefore, the asymptotic expressions for the Bessel functions with large arguments become now relevant. For large z,|z|  1, Bessel

023007-6

STABILITY ANALYSIS OF A THIN LIQUID FILM ON . . .

PHYSICAL REVIEW E 90, 023007 (2014)

functions of nth order satisfy the following relationships:

√ −z πe ez In (z) ∝ √ , Kn (z) ∝ √ , 2π z 2z

and the evolution equation (31a) reduces to  2  √ √ 3B  hx [JL ( 2h) + DL ( 2h)] + h3 hx + h3 hxxx = 0. ht + C x

(33a)

Here JL (h) and DL (h) denote the corresponding limits of J (n) and D(n) functions of Eq. (31a) in the limit a → 0, [6 sin h sinh h − 3h(cos h sinh h + sin h cosh h) − 8 sin (h/2) sinh (h/2)] , 4(cos h + cosh h)2 2 sin (h/2) sinh (h/2) . DL (h) = (cos h + cosh h)2 JL (h) =

Note that in this limiting process, the overall contribution of the remaining terms of L1 is zero. Also note that in the limiting process of a → 0 for the function J (r,n), we used a transformation r → 1 + az. Equation (33) is identical to the equation derived by Shklyaev, Alabuzhev, and Khenner [20], and further investigated by Benilov and Chugunova [21]. Equation (31a) contains three independent dimensionless parameters, namely a, ϒ(=B 2 /C), and . We note that if the properties of the liquid and the average film thickness are fixed, thus C remains fixed, and the behavior of the system is investigated in the context of the three parameters a, B, and .

(33b)

analysis derived by Goren [9] in the case of a negligible inertia in terms of Bessel functions without an assumption of a long-wave limit. To enable matching between the values, as given in Eqs. (17)–(19) there, the following substitutions have to be made: S → C/kc , ka → k/kc , ks → k/a, and N → λ/kc2 . Figure 1 displays the variation of the normalized growth with the normalized wave number kn = kkc , rate λn = a(1+a)λ Y (n0 ) where kc is taken from Eq. (37). The growth rate arising from our theory and given by Eq. (36) is shown by the solid

V. LINEAR STABILITY ANALYSIS

In this section, we carry out the linear stability analysis of the system based on the evolution equation (31a). The base state considered is n = n0 ≡ 1 + a corresponding to the undisturbed cylindrical interface of a uniform thickness h¯ = 1. Equation (31a) is linearized around the base state via n = n0 + εy,ε  1 to yield yt + yζ ζ L1 (n0 )/n0 + yζ ζ ζ ζ L2 (n0 )/n0 = 0.

[L1 (n0 ) − k L2 (n0 )] . n0

0.8

(34)

Substituting a normal-mode disturbance y∼ exp[ikζ + λ(k)t] with the wave number k and growth rate λ, both dimensionless, results in the dispersion relation λ(k) = k 2

1

0.6

λn

0.4

2

(35) 0.2

In the case of an unforced system, B = 0, Eq. (35) may be rewritten as   2 Y (n0 )k 2 a 2 (36) λ(k) = −k , a(1 + a) 1 + a and one can infer that the film is unstable for k < kc and stable for k > kc , where the cutoff wave number is determined [16] by a kc = . (37) 1+a To validate our treatment of the capillary term discussed above we wish to begin with the comparison between the dispersion relation derived above and given by Eq. (36) in the particular case of the unforced system, B = 0, and the dispersion relation arising from a general linear stability

0

0

0.2

0.4

0.6

kn

0.8

1

FIG. 1. (Color online) Variation of the normalized dimensionless λ , with the normalized wave number kn = kkc , growth rate λn = λmax as given by Eq. (35) and Eqs. (17)–(19) of Goren [9] in the static case, B = 0. The solid curve represents the long-wave result obtained from Eq. (36), whereas the dotted, dashed, and dot-dashed curves, respectively, correspond to the general results of Goren [9] for a = 0.1, 0.2, and a = 0.3. As a result of normalization for the wave number and the growth rate, the curves for the two different liquids considered here collapse on a single respective curve depending only on a.

023007-7

SELIN DURUK AND ALEXANDER ORON

PHYSICAL REVIEW E 90, 023007 (2014) 0.1

5 0.1

0.08

4 2

0.05

0.06

3 2

S 0

S

1

0.04

1 2

3

0.02

3 -0.05

4 0

4

5 6

5 0

5

10

15

20

Ω

25

30

35

-0.02 0.1

40

0.2

0.3

a

0.4

0.5

FIG. 2. (Color online) Left panel: variation of S(n0 ,) with  for various values of a = 0.1, 0.2, 0.3, 0.4, and 0.5, as represented by curves 1–5, respectively. Right panel: Variation of S(n0 ,) with a for various values of  = 1, 4, 8, 9, 9.5, and 10, as represented by curves 1–6, respectively.

curve. The corresponding curves obtained from Goren [9] are presented by broken curves for varying a. A very good agreement between the two theories is clearly observed for small a. For instance, the difference between the two theories for a = 0.3 is less than 10%. The difference between the two theories increases with an increase in a around the wave number corresponding to the fastest growing linear mode. It is important to mention that the curves shown in Fig. 1 are independent of liquid properties in the given scaling Eq. (30). In what follows, we carry out the linear stability analysis of the forced system. In order to reduce a number of the dimensionless parameters of the problem, we focus on two specific liquids. These are water-glycerin solutions whose material properties are as follows [30]: (i) 33% glycerin by weight, with ρ = 1075 kg/m3 ; ν = 2.9 × 10−5 m2 /s; σ = 69 × 10−3 kg/s2 ; in normal gravity the capillary length is c ≈ 2.5 mm; (ii) 80% glycerin by weight, with ρ = 1210 kg/m3 ; ν = 6.9 × 10−5 m2 /s; σ = 66 × 10−3 kg/s2 ; in normal gravity the capillary length is c ≈ 2.3 mm. In the beginning of Sec. II, we imposed a condition in terms of the Bond number for validity of neglecting the gravity effects. Rewriting the inequality Bo < 0.5 in terms of the parameters a, h0 , and c for the data used here, results in a > 0.1. On the other hand, we are interested in relatively thin films in comparison with the substrate radius, a < 0.5. Hence the range for a considered here will be 0.1  a  0.5. Once the external forcing is present, B = 0, the behavior of L1 (n0 ) is determined by the geometric parameter a and the forcing parameters ϒ and . As seen in Fig. 2, the function S(n0 ,) = J (n0 ,) + a 2 D(n0 ,), which represents a forcing-related component of the function L1 (n0 ) may change its sign as a function of  for fixed values of a. The function S(n0 ,) depends on two independent variables, namely n0 and ; thus it is instructive to look at its variation in two different cross sections. As displayed in the left panel of Fig. 2, S(n0 ,) increases for small fixed a until it reaches

its maximal value at around   2.9, and then decreases to its root around   9.45, remaining thereafter negative in the range investigated here,   40. A similar feature is also observed in the right panel of Fig. 2, as S(n0 ,) first increases remaining positive for all a for smaller values of , and decreases thereafter becoming negative for larger values of ,  > 9.45, for all a. However, to qualitatively change the stability properties of the system by changing the sign of L1 (n0 ), the value of S(n0 ) must exceed a certain positive threshold in order to overcome the effect of the curvature contribution to the system instability expressed by L1 (n0 ). The effect of vibration on the film evolution is determined by the ϒ term in Eq. (31a) which consists of three components. These are the two elements of the inertial terms arising from the pulsating flow 

uu x  and 

w u r , and the term arising from the coupling between the axial component of the pulsating flow field

u and the pulsating component of the interface location

h. We will refer below to these as terms 1, 2, and 3, respectively. They are each related, respectively, to the f1 , f2 contributions to the J term, and the D term in Eq. (29a). The sum of terms 1 and 2 yields a combined vibrationally induced correction to the axial velocity component u¯ of the averaged flow field, Eq. (27), with respect to the pressure-gradient term in the case of the unforced flow subject to the capillary instability. The sum of all the three terms represents the function S(n0 ,). We find that in a typical situation the contribution of term 3 is small for small a with respect to the others, associated with the fact that it contains a 2 as a factor. We also find that the contributions of terms 1 and 3 to the function S are typically positive, i.e., destabilizing, whereas that of term 2 is negative and thus stabilizing. Therefore, we conjecture that the major role in a possible stabilization of the axisymmetric flow is provided by the radial component of the pulsating flow field w

. Its presence weakens and may also revert the flow induced by the destabilizing axial pressure gradient due to capillarity.

023007-8

STABILITY ANALYSIS OF A THIN LIQUID FILM ON . . .

PHYSICAL REVIEW E 90, 023007 (2014)

0.2

0.1

1

3 0

2 -0.1

-0.2 0.1

0.2

0.3

0.4

a

0.5

FIG. 3. (Color online) Variation of the three separate components contributing to the film evolution by vibration. Curves 1, 2, and 3 correspond, respectively, to the f1 and f2 components in the J (n) term, and the D(n) term in Eq. (29a). The combination of these terms represents the vibrational component of the function L1 (n) in Eq. (31a). The solid and dashed curves correspond to  = 8 and  = 9.6, respectively.

Figure 3 presents the variations of the three aforementioned components of the function S(,n0 ) with the parameter a for 0.4

0.3

0.006

7

1.2 0.8

two cases of  = 8 and  = 9.6 corresponding to S > 0 and S < 0, respectively. All of the three terms, 1, 2, and 3, exhibit a monotonic behavior with respect to a. As for their variation with , their absolute values decrease with an increase in  for  > 4. We reiterate that the instability of the system is fed by the combination of the capillary instability and the effect of vibration via the function S discussed above. In order to neutralize the capillary instability one needs to manipulate the vibration parameters so that the function S becomes negative with a sufficiently large absolute value. It thus follows that, if S < 0, an increase in the parameter ϒ which is proportional to the squared dimensionless forcing amplitude B enhances the tendency for stabilization. Relative variation of the terms 1, 2, and 3 with a enables the total instability term L1 containing their sum with the capillary contribution related to the azimuthal component of the interfacial curvature to produce a dip in the curve of S as a function of a possibly protruding into the negative domain with an increase in . There, stabilization of the axisymmetric flow takes place. This dip typically widens with an increase of B (via an increase in ϒ). Figure 4 presents the results of linear stability analysis in the case of a 33% water-glycerin solution film with h0 = 0.2 mm, a fixed dimensionless forcing amplitude of B = 5 (b = 1 mm), and for various values of dimensionless forcing frequency in the range of 1.3793    9.7931 corresponding to the frequencies 1000 Hz  ω  7100 Hz. In this case, the value of ϒ is moderate, ϒ = 1.6378. The left panel of√ Fig. 4 presents the variation of the cutoff wave number kc = L1 (n0 )/L2 (n0 ) [where L1 (n0 )  0] and kc = 0, where L1 (n0 )  0, with a for the chosen parameter sets. The dashed curve there represents the cutoff wave

6

4×10

8

0.4

-5

7 86

5

0

0

0.5

0.2

5

0.003

0 0.4 4

kc 0.2

λ

4

3

0 2 1

3

0.1

1

2

-4×10-5

0 0.1

0.2

0.3

a

0.4

0.5

0

0.1

k

0.2

FIG. 4. (Color online) Linear stability analysis for a 33% water-glycerin solution film with a thickness of h0 = 0.2 mm. The dashed curves correspond to the static case b = 0. Left panel: variation of the cutoff wave number kc with the parameter a for a fixed forcing amplitude b = 1 mm (B = 5). In this case ϒ = 1.6378. Solid curves 1–8 show, respectively, the cases of ω = 7100, 7000, 6950, 6900, 6800, 6000, 2000, and 1000 Hz ( = 9.7931, 9.6552, 9.5862, 9.5172, 9.3793, 8.2785, 2.7586, and 1.3793). Right panel: the corresponding growth rate λ as a function of the wave number k for a = 0.3 and the values of b and ω as presented for the left panel. 023007-9

SELIN DURUK AND ALEXANDER ORON

PHYSICAL REVIEW E 90, 023007 (2014)

5

0.3

2×10-6 4 3

2

0.2

kc

λ

1

5 0.1

4

3

2

0

2 1

1

6 -2×10-6

0 0.1

0.2

0.3

a

0.4

0

0.5

0.05

0.1

k

0.15

0.2

FIG. 5. (Color online) Linear stability analysis for a film of 33% water-glycerin solution with a thickness of h0 = 0.2 mm. The dashed curves correspond to the static case, b = 0. Left panel: variation of the cutoff wave number kc with the parameter a for two values of the fixed forcing frequency. Dot-dashed curves 1 and 2 correspond to ω = 6000 Hz with b = 0.1 mm and 0.2 mm, respectively. Solid curves 1–6 correspond, respectively, to ω = 7000 Hz, with the forcing amplitudes b = 1, 0.9, 0.7, 0.5, 0.4, and 0.1 mm (B = 5, 4.5, 3.5, 2.5, 2, and 0.5,  = 9.6552). Right panel: the corresponding growth rate λ(k) for a = 0.2. Curves 1–5 display, respectively, the cases of ω = 7000 Hz, b = 0.9, 0.7, 0.5, 0.4, and 0.1 mm (B = 4.5, 3.5, 2.5, 2, and 0.5,  = 9.6552).

number kc for the case of the static cylinder, B = 0, as given by Eq. (37). With an increase of the forcing frequency ω with a fixed forcing amplitude b, the instability domain grows as indicated by an increase in kc ; see curve 5 (ω = 6800 Hz) in the left panel of Fig. 4 and curves 7 and 8 (ω = 2000 Hz and ω = 1000 Hz) in the inset there. With a further increase in , the range of unstable modes begins to shrink; see curve 6 (ω = 6000 Hz) in the inset. More increase in the forcing frequency ω yields stabilization of the forced axisymmetric flow with respect to the static one, as shown by curves 4–1 (ω = 6900 Hz to ω = 7100 Hz) in the left panel of Fig. 4. Curves 1–3 (ω = 7100 Hz to ω = 6950 Hz) in the left panel of Fig. 4 branch off zero at the value of a = a0 depending on , so that the axisymmetric flow is stable for 0.1 < a < a0 (except for the case corresponding to curve 5 in Fig. 4, where instability takes place for all a) and unstable for a > a0 . The value of a0 increases with ; thus an increase in the forcing frequency stabilizes the axisymmetric flow for a wider range of a. The cutoff wave number kc , when positive, increases with a for a fixed forcing frequency. When the system is unstable under forcing, the range of linearly unstable modes tends to be narrower in the unforced case for larger a. Also, for larger a the range of linearly unstable modes shrinks with an increase in  when a is fixed. Note that for ω > 7200 Hz ( > 9.9310) the film is stable for all a in the considered range. The right panel of Fig. 4 shows the growth rate λ as a function of the wave number k for a fixed value of the forcing parameter b = 1 mm, a fixed value of the geometrical parameter a = 0.3, and for the values of the forcing frequency ω used in the left panel of Fig. 4. Again, the dashed curve

corresponds to the case of a static cylinder, B = 0, with the same value of a. When the system is unstable, the growth rate of the disturbances tends to increase when forcing is applied for lower frequencies and to decrease for higher ones. It is clearly seen in the right panel of Fig. 4 that destabilization and stabilization of the axisymmetric flow is accompanied, respectively, by an increase and decrease in the growth rate for every value of the disturbance wave number k. Figure 5 presents in its left panel the variation of the cutoff wave number kc with the parameter a for two fixed forcing amplitudes and for various forcing frequencies. As in the previous case, the dashed curve shows the cutoff wave number for the static case, B = 0. In the case of a lower forcing frequency, ω = 6000 Hz, forcing with amplitudes up to the mean thickness of the film results in a system destabilization in terms of the width of the instability range. In the case of a larger forcing frequency, ω = 7000 Hz, forcing with small amplitudes up to a film thickness practically does not affect the cutoff wave number kc relative to the static case. A further increase in the forcing amplitude leads to the emergence of a stability window similar to the previous case addressed in Fig. 4. The stability window widens in terms of a and the instability domain narrows with an increase in the forcing amplitude ω; see the left panel of Fig. 5. The growth rate λ presented in the right panel of Fig. 5 for a = 0.2 decreases for all k, as the amplitude b increases. In the case of a more viscous liquid, an 80% water-glycerin solution with a much larger value of ϒ, the topology of the relevant curves is different. Figure 6 displays the variation of the cutoff wave number kc with a for varying frequency

023007-10

STABILITY ANALYSIS OF A THIN LIQUID FILM ON . . .

PHYSICAL REVIEW E 90, 023007 (2014)

2×10-5

0.3

0.2

kc

λ

0

1 2

0.1

1

3

2

3

4

3

4

4 -5

0 0.1

0.2

0.3

0.4

a

0.5

-2×10

0

0.1

0.2

k

FIG. 6. (Color online) Linear stability analysis for a film of 80% water-glycerin solution with a thickness of h0 = 0.2 mm. In this case, ϒ = 10.9106. The dashed curves correspond to the case of a static system b = 0. Solid curves 1–4 show, respectively, the cases of ω = 16380, 16385, 16390, and 16395 Hz ( = 9.4957, 9.4986, 9.5014, and 9.5043). Left panel: variation of the cutoff wave number kc with the parameter a for a fixed forcing amplitude b = 1 mm (B = 5). Right panel: variation of the growth rate λ(k) with the disturbance wave number for a fixed a = 0.3.

of forcing with a fixed forcing amplitude. For small forcing frequencies, the system is destabilized beyond the instability of the static case (not shown in Fig. 6). However, as the forcing frequency becomes sufficiently large, a dip emerges in the variation of kc with a which expresses destabilization of the system (broadening of the instability range and increase of the growth rate, as seen in the right panel of Fig. 6) with respect

to the static one for smaller a and its stabilization (shrinking of the instability domain and decrease of the growth rate) with respect to the static case for larger but still small a. For even higher forcing frequencies a window of a full stabilization of the axisymmetric flow emerges in the intermediate range of a. This stability window widens in both directions as the forcing frequency ω increases.

2×10-5

0.3

1

2

0.2

1

2

kc

λ

0 3 4

0.1 3

-2×10

-5

4 6 0 0.1

5

5 0.2

0.3

a

6 0.4

0.5

0

0.1

k

0.2

FIG. 7. (Color online) Linear stability analysis for a film of a 80% water-glycerin solution with a thickness of h0 = 0.2 mm. The dashed curves correspond to the case of a static system b = 0. Left panel: variation of the cutoff wave number kc with the parameter a for a fixed forcing frequency ω = 16390 Hz ( = 9.5014). Solid curves 1–6 show, respectively, the effect of varying forcing amplitude for b = 0.1, 0.6, 0.9, 0.95, 0.98, and 1 mm (B = 0.5, 3, 4.5, 4.75, 4.9, 5). Right panel: variation of the growth rate λ(k) with the disturbance wave number for a fixed a = 0.3. Curves 1–4 show, respectively, b = 0.1, 0.6, 0.95, and 1 mm (B = 0.5, 3, 4.75, and 5). 023007-11

SELIN DURUK AND ALEXANDER ORON

PHYSICAL REVIEW E 90, 023007 (2014)

Figure 7 presents a variation of the cutoff wave number kc and the growth rate λ for a fixed forcing frequency with a varying forcing amplitude. The forcing leads to stabilization of the axisymmetric flow (both shrinkage of the instability range and a decrease in the growth rate) with an increase in the amplitude b. A full stabilization window develops in the intermediate range of a. The full stability window becomes wider in terms of a in both directions with an increase in the forcing amplitude. Finally, we emphasize that a difference in the topology of the neutral curves shown in Figs. 4 and 5 on one hand and Figs. 6 and 7 on the other hand which exhibit, respectively, stabilization windows on the left and in the middle of the displayed range of a is only apparent. It arises because of the constraint on a equivalent to the requirement of a small Bond number discussed in the beginning of this section. Were the results of Figs. 4 and 5 formally extended into the range of smaller a, they would exhibit the same topology as in Figs. 6 and 7, namely, with a stability window in the middle of the relevant a domain. VI. BIFURCATION ANALYSIS

In this section, we carry out the weakly nonlinear analysis based on Eq. (31a). Linear stability analysis presented above yields the critical value of ϒ = ϒc , ϒc =

  a 2 Y (n0 ) k 2 − 1+a a[J (n0 ) + a 2 D(n0 )]

.

(38)

The following expansions are now introduced into Eq. (31a): T = δ 2 t, n = n0 + δn1 + δ 2 n2 + δ 3 n3 + · · · ,

with C representing an unknown complex amplitude function. At third order we obtain n0 nT + L2 (n0 )n3,ζ ζ ζ ζ + L1 (n0 )n3,ζ ζ + sn1,ζ ζ S(n0 ) + v1 (n0 )(n1 n2 )ζ ζ + L2 (n0 )(n1,ζ ζ ζ n2 + n2,ζ ζ ζ n1 )ζ + 12 (n21 n1,ζ )ζ L1 (n0 ) + 12 (n21 n1,ζ ζ ζ )ζ L2 (n0 ) = 0. (44) Eliminating the secular term of the equation at third order in δ yields the following Landau equation in terms of the amplitude function A: (1 + a)Aτ = sS(n0 )A − |A|2 A, where =

VII. SUMMARY

(40)

in the domain of 0  ζ  l/ h0 = −1 . Looking for periodic solutions with a zero mass in order to not violate the mass conservation with respect to the base state, we find [since L1 (n0 ) = k 2 L2 (n0 )] (41)

where k = 2π is the fundamental wave number and A(T ) is a complex amplitude function to be determined below. At second order we obtain L2 (n0 )n2,ζ ζ ζ ζ + L1 (n0 )n2,ζ ζ

n2 = C(T )eikζ +

A(T )2 (L1 − L2 ) 2ikζ e + c.c., 6v1

(45b)

(39)

where δ  1 is a small deviation from criticality and s = ±1. Substituting Eqs. (39) into Eq. (31a) results in a hierarchy of problems at each order of δ. At first order, we obtain a fourth-order linear ordinary differential equation,

= −[L1 (n0 )(n1 n1,ζ )ζ + L2 (n0 )(n1 n1,ζ ζ ζ )ζ ],

L2 (n0 ) − L1 (n0 ) 2

is the Landau constant whose sign determines the bifurcation type. We find that for all regimes considered here,  is always negative and therefore the bifurcation is subcritical. Since near the threshold of linear instability, S(n0 ) < 0, the instability takes place for negative s, i.e., with a decrease of ϒ from its critical value. Subcriticality of the bifurcation in the given system implies that near the onset of the linear instability, the system is expected to propel toward rupture when the system is unstable. However, farther away from criticality, a full numerical investigation of Eq. (31a) is necessary to be carried out in order to investigate the nonlinear dynamics of the forced system.

ϒ = ϒc + sδ ,

n1 = A(T )eikζ + c.c.,

[7L2 (n0 ) − L1 (n0 )][L1 (n0 ) − L2 (n0 )] 6L1 (n0 ) +

2

L2 (n0 )n1,ζ ζ ζ ζ + L1 (n0 )n1,ζ ζ = 0,

(45a)

(42) (43)

In this paper we have investigated stability of an axisymmetric thin liquid film on a horizontal circular cylinder undergoing fast axial harmonic oscillation. Using the methods of long-wave theory combined with the multiscale expansions, separation of the relevant fields into fast and slow components, and averaging over the fast time, we have derived a nonlinear evolution equation describing the spatiotemporal dynamics of the average film thickness in a slow time scale. Based on the assumption of a small film parameter and scalings used here, we reckon that the regime of a high-frequency limit is achieved already if the condition   2 is satisfied. In the limiting case of zero forcing amplitude, our evolution equation reduces to Eq. (25) derived in Weidner, Schwartz, and Eres [14], and later in Haimovich and Oron [16]. It is further reduced in the limit of a small parameter a to the equation derived by Hammond [10]. In the case of a forced system, B = 0, our derived evolution equation yields in the limit of a large cylinder radius the equation derived for a forced film on a planar substrate [20,21]. Linear stability analysis of the system has been carried out in the framework of the derived evolution equation. The results obtained from the linear stability analysis have been compared with the behavior of the unforced system which is intrinsically unstable, in terms of the cutoff wave number and

023007-12

STABILITY ANALYSIS OF A THIN LIQUID FILM ON . . .

PHYSICAL REVIEW E 90, 023007 (2014)

the growth rate of disturbances. For specific fluids chosen here as water-glycerin solutions and for a fixed mean film thickness, we have found that it is possible to stabilize an axisymmetric film coating on a cylindrical surface in a certain range of a for normal modes of all wave numbers. This range of stabilization depends on a set of forcing parameters. One of the reasons for such a stabilization is certainly due to a liquid inertia associated with the forcing imparted on the film via the boundary. On the other hand, the same parameter set could be also destabilizing for a cylinder of a different radius. Hence fast harmonic axial forcing of a cylinder may be a feasible tool for manipulating an unstable character of a liquid film on a static cylinder, making it either stable or even more unstable by increasing the growth rate of disturbances and/or by widening the range of linearly unstable modes. It is also important to note that the results of the linear stability analysis carried out here are valid for small wave numbers k of the normal modes with respect to which stability properties of the system are studied.

We have also carried out a weakly nonlinear analysis of the system and investigated the character of the bifurcation from the equilibrium by monitoring the sign of the relevant Landau constant . The bifurcation has been found subcritical. This implies that beyond and near the stability threshold the film is expected to rupture. Finally, a thorough time-dependent numerical investigation of Eq. (31a) is not attempted in this paper. Its results will be reported elsewhere to support our conclusions made based on both linear stability and bifurcation analyses, and to explore the nonlinear dynamics of the system farther away from its stability threshold.

[1] A. Stephenson, Manchester Memoirs 52, 1 (1908). [2] P. L. Kapitza, Zh. Eksp. Teor. Fiz. 21, 588 (1951) (in Russian). Also, translated in P. L. Kapitza in Collected Papers, edited by T. den Haar (Pergamon, London, 1965), Vol. 2, p. 714. [3] P. Brunet, J. Eggers, and R. D. Deegan, Phys. Rev. Lett. 99, 144501 (2007). [4] L. Moldavsky, M. Fichman, and A. Oron, Phys. Rev. E 76, 045301 (2007). [5] X. Noblin, R. Kofman, and F. Celestini, Phys. Rev. Lett. 102, 194504 (2009). [6] C.-S. Yih, J. Fluid Mech. 31, 737 (1968). [7] A. C. Or, J. Fluid Mech. 335, 213 (1997). [8] D. Halpern and A. L. Frenkel, J. Fluid Mech. 446, 67 (2001). [9] S. Goren, J. Fluid Mech. 12, 309 (1962). [10] P. S. Hammond, J. Fluid Mech. 137, 363 (1983). [11] P. A. Gauglitz and C. J. Radke, Chem. Eng. Sci. 43, 1457 (1988). [12] A. L. Yarin, A. Oron, and P. Rosenau, Phys. Fluids 5, 91 (1993). [13] M. Johnson, R. D. Kamm, L. W. Ho, A. Shapiro, and T. J. Pedley, J. Fluid Mech. 233, 141 (1991). [14] D. E. Weidner, L. W. Schwartz, and M. H. Eres, J. Colloid. Interface Sci. 187, 243 (1997).

[15] [16] [17] [18] [19]

ACKNOWLEDGMENTS

The research was partially supported by the European Union via the FP7 Marie Curie scheme [PITN-GA-2008214919 (MULTIFLOW)]. The authors acknowledge many fruitful discussions with Dr. Sergey Shklyaev.

[20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

023007-13

B. Reisfeld and S. G. Bankoff, J. Fluid Mech. 236, 167 (1992). O. Haimovich and A. Oron, Phys. Fluids 22, 032101 (2010). O. Haimovich and A. Oron, Phys. Rev. E 84, 061605 (2011). O. Haimovich and A. Oron, Phys. Rev. E 87, 052403 (2013). S. Shklyaev, M. Khenner, and A. A. Alabuzhev, Phys. Rev. E 77, 036320 (2008). S. Shklyaev, A. A. Alabuzhev, and M. Khenner, Phys. Rev. E 79, 051603 (2009). E. S. Benilov and M. Chugunova, Phys. Rev. E 81, 036302 (2010). V. Lapuerta, F. J. Mancebo, and J. M. Vega, Phys. Rev. E 64, 016318 (2001). U. Thiele, J. M. Vega, and E. Knobloch, J. Fluid Mech. 546, 61 (2006). A. Oron, S. H. Davis, and S. G. Bankoff, Rev. Mod. Phys. 69, 931 (1997). R. V. Craster and O. K. Matar, Rev. Mod. Phys. 81, 1131 (2009). J. R. de Bruyn, Phys. Fluids 9, 1599 (1997). R. V. Craster and O. K. Matar, J. Fluid Mech. 553, 85 (2006). E. Novbari and A. Oron, Phys. Fluids 21, 062107 (2009). E. Novbari and A. Oron, Phys. Fluids 23, 012105 (2011). M. F. G. Johnson, R. A. Schluter, M. J. Miksis, and S. G. Bankoff, J. Fluid Mech. 394, 339 (1999).