Thin film hydrodynamics of second- and third-order fluids confined in a Hele–Shaw cell

International Journal of Engineering Science 43 (2005) 1388–1408 www.elsevier.com/locate/ijengsci

Thin film hydrodynamics of second- and third-order fluids confined in a Hele–Shaw cell Fethi Kamısßlı

*

Department of Chemical Engineering, Faculty of Engineering, Firat University, 23119 Elazig, Turkey Received 25 October 2004; received in revised form 28 April 2005; accepted 23 May 2005

Abstract The thin film coating developed by penetration of a long gas bubble through a Hele–Shaw cell initially filled with a viscoelastic fluid is theoretically investigated using the third-order fluid model. The work presented in this paper is directed at identifying effects of fluid elasticity on the liquid fraction deposited on the walls of a Hele–Shaw cell. The singular perturbation method is used to determine the residual liquid film thickness of a viscoelastic fluid on the wall of a planar geometry when displaced by another immicible fluid. Inner and outer expansions are developed in terms of small parameters Ca1/3, d (d = De/Ca1/3) and d2. The method of matched asymptotic expansions is used to match the inner and outer solutions by means of transition region between the advancing meniscus and the entrained film where the fluid rheology has its greatest effect. An understanding of the role of fluid elasticity comes from a consideration of both the second- and thirdorder fluid models, a detailed analysis indicates that the residual liquid film thickness of the viscoelastic fluid tends to decrease while the pressure drop across the bubble front tends to increase as the fluid becomes more viscoelastic. While the force due to the present of a normal stress gradient in the flow direction was evaluated from the second-order fluid model and acts to decrease the film thickness and increase the pressure jump, whereas the force due to shear stress thinning was evaluated from the third-order fluid model and tends to increase the film thickness and decrease the pressure jump. Theoretical results presented in this paper are in agreement with other the experimental and theoretical studies. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: Second- and third-order fluids; Two-phase flow; Perturbation method

1. Introduction Gas-assisted displacement processes are commonly encountered in a wide range of industrial applications such as in enhanced oil recovery, production of hollow fiber membranes, coating of ceramic monoliths for the manufacture of catalytic converters, gas–liquid contacting (i.e., bubble column, vapor–liquid contactors or absorbers, reboilers, spraying system, etc.), activated sludge processes and in polymer processes such as gas-assisted injection molding. *

Tel.: +90 424 23700; fax: +90 424 5526. E-mail address: fkamisli@firat.edu.tr

0020-7225/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijengsci.2005.05.011

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

1389

Many engineering applications involve the entrainment of liquid films on solid objects and prediction of film thickness is of interest in pickling, washing, lubrication and especially coating operations. In gas-assisted injection molding a mold is partially filled with a polymer melt. This is followed by injection of pressurized gas through the polymer melt to fill the mold cavity. Therefore, hollow parts can be produced with higher strength to weight ratio, improved surface finish and lower molded in stresses. When a less viscous fluid displaces a more viscous one from the gap between two closely spaced parallel plates, the interface develops a tongue-like shape with the less viscous fluid penetrating into the more viscous one. Similarly, when air is forced into one end of a circular tube containing a viscous liquid, it forms a roundended column or bullet-like shape which travels down the tube forcing some liquid out at the far end and leaving a fraction of the liquid, m, in the form of an annular layer covering the wall. Gas-assisted displacement has been studied for more than seventy years. One of the earlier experimental studies was performed by Fairbrother and Stubbs [1]. The experimental results obtained from penetration of an air bubble through a Newtonian fluid confined in a circular tube were formulated with an empirical correlation between fraction of the liquid and capillary number as follows:
gu 1=3 ub
u b m¼ ¼ 1:0Ca1=3 ¼ ; ð1Þ ub r where g is the viscosity of the liquid phase, r the surface tension at the interface, u the displaced fluid velocity and ub the bubble velocity. This correlation was found to be satisfactory for capillary numbers between 10
3 and 10
2. The fractional coverage, m, is defined as the fraction of the tube cross-sectional area coated with fluid after a long bubble penetration, and the capillary number is a dimensionless group which represents the ratio of viscous force to surface tension in the gas-assisted displacement process. The limit of large surface tension corresponds to the case of the capillary number approaching zero. The above experimental study was extended by some investigators [2–4] to higher viscosity of liquids (and thus, capillary number) to determine the fractional coverage of liquid deposited on the wall of the tube as a function of capillary number. These studies indicated that the fractional coverage of liquid asymptotically increased to the values of 0.56–0.60. A theoretical analysis of the gas-assisted displacement of a Newtonian fluid confined in either a circular tube or a rectangular channel (Hele–Shaw cell) was undertaken by many investigators [5–9] using the method of matched asymptotic expansions. The idea behind this theoretical treatment is that for sufficiently small Ca the viscous stresses appreciably modify the static profile of the bubble only very near to the wall. In this region the lubrication approximation gives a good description of the flow field and of the interface profile. In the center of the capillary, the static profile is valid and there is a region of overlap in which the two solutions are matched. Using the lubrication approximation which requires quasi-unidirectional flow in the thin liquid film and assuming the slope of the fluid–fluid interface to be small, it can be shown that the velocity profile is parabolic. The boundary conditions for steady flow are the no slip condition at the capillary wall, and tangential stress equal to zero at the fluid interface. The bubble is assumed to be inviscid resulting in a constant pressure within the bubble. The pressure in the liquid film is given by the pressure drop across the interface which is approximated by the Young–Laplace equation. The isothermal gas-assisted displacement of a Newtonian liquid in a channel of square cross-section was studied by Kolb and Cerro [10] to determine the fractional coverage of liquid that was found to be 0.64 at high capillary numbers. They observed that the hollow core developed by the less viscous fluid penetrating the more viscous one takes either the shape of the square cross-section of the channel or a bullet-like shape depending on the value of the capillary number. If the capillary number is less than 0.1, the long bubble (hollow core) conforms to the shape of the square cross-section and thus, less liquid deposits on the walls. Otherwise (Ca > 0.1), it forms a circular hollow core and thus, the more liquid deposits on the walls. Furthermore, in recent years there have been several studies [11–14] on the gas-assisted displacement of non-Newtonian fluids, not only because of their technological significance but also in the interesting mathematical features presented by equations governing the flow. Huzyak and Koelling [12] experimentally investigated the effect of liquid elasticity on the residual liquid film thickness of viscoelastic fluids. The effect of

1390

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

fluid elasticity on the liquid film thickness of a viscoelastic fluid confined in a Hele–Shaw cell was theoretically examined using an Oldroyd-B constitutive equation [11]. The planar geometry or Hele–Shaw cell consists of two closely separated parallel plates having a distance of 2d between them. The sides of this rectangular channel are at a distance of 2Z0 apart where d  Z0, parameter k is defined as (thickness of gas bubble)/(distance between the plates). In a rectangular channel the thickness of the tongue-like shape is 2 kd and its width is 2 kwZ0, where the parameter kw is equal to (width of the bubble)/ (width of the rectangular channel). The determination of the value of k and kw in a planar geometry has been a subject of much interest. The gas-assisted displacement of a Newtonian fluid in a Hele–Shaw cell has been experimentally studied by many authors [15–18] to determine kw as a function of capillary number, Ca, for different cell aspect ratios, Z0/d. The elasticity of the fluid seems to give rise to interesting phenomena when a viscoelastic fluid is displaced by an inert gas [11,12]. Bauman et al. [19] have been also observed experimentally that the addition of a dilute polymer has a destabilizing effect on two-roll coating flow, giving rise to lower critical values above which ribbing instability occurs. The stability of the interface for two-phase immiscible displacement of Newtonian fluids in a planar geometry was also investigated theoretically as well as experimentally [15]. The effect of fluid elasticity on the residual liquid film thickness was examined experimentally by Huzyak and Koelling [12]. They observed that at the low Deborah number (De < 1), the value of reduced fractional coverage, mR = m/mN, (m: fractional coverage of viscoelastic fluid, mN: fractional coverage of Newtonian fluid) is between 0.95 and 1.05. The value of reduced fractional coverage less than unity indicates that the film thickness of viscoelastic fluid is less than that of a Newtonian one. A value of reduced fractional coverage mR = 1.0 implies that the fractional coverage of viscoelastic fluid is identical with that of the Newtonian fluid at an equivalent capillary number. The fractional coverage of viscoelastic fluids continues to increase with Deborah number for all De P 1 and approaches to a value that is 30% greater than the Newtonian result at high Deborah numbers (De ffi 5). The theoretical analysis [11] also indicated that the film thickness of the viscoelastic fluid described by an Oldroyd-B equation tends to decrease and the pressure drop at the interface tends to increase as the fluid becomes more viscoelastic. Although the secondorder fluid model is able to predict the normal stress differences, which are characteristic of non-Newtonian fluids, it does not take into account the shear thinning and thickening phenomena that many real fluids exhibit. The second-order fluid is the second-order asymptotic approximation about the state of rest of a viscoelastic fluid. Hence, the applicability of the second-order fluid is restricted to slow and slowly varying flow fields. The third-order fluid model represents a further, although inconclusive, attempt towards a comprehensive description of the properties of viscoelastic fluids. With the above experimental and theoretical observations in mind, the elasticity of the fluid seems to give rise to interesting phenomena when a viscoelastic fluid is displaced by an inert gas bubble. Therefore, to observe the differences among different constitutive equations for viscoelastic fluids, the model in the present paper is chosen as a third-order fluid. The task here is to find an exact analytical solution for the gas-assisted displacement of third-order fluid confined in a Hele–Shaw cell. The fluid considered is displaced by an inert gas in a thin gap which consists of two parallel plates separated by a small distance of 2d. It is the purpose of this paper to examine the effect of fluid elasticity on the thin film creeping flow developed by a less viscous fluid penetrating a more viscous one in a Hele–Shaw cell initially filled with a viscoelastic fluid that is described by the third-order fluid model. Although the viscoelastic fluid is modeled by the third-order fluid model, the resulting governing differential equation contains some variables that are to be determined from the first-, second- and third-order fluid models as will be seen later. Thus, the first-, secondand third-order fluid models are used together to evaluate the effect of fluid elasticity on the residual thickness of the third-order fluid and the pressure drop during the gas-assisted displacement process. In the lubrication approximation, it has been assumed that the bubble moves with sufficiently small velocity, and the interface is almost parallel to the wall of channel. Thus, the slope of the fluid–fluid interface is small and the unidirectional flow takes place in the thin liquid film. For an analytical solution, it was assumed that the residence time in the transition region is much longer than the relaxation time of polymer. In that case, the polymeric fluid behaves essentially as a Newtonian fluid. The effect of a finite relaxation time which is long enough to disturb the local

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

1391

flow and stress fields away from their Newtonian values can be seen by examining the dynamics at O(d = De/Ca1/3). Since it is assumed that the thickness of the gap between the two plates is very small, it is possible to derive two dimensional field equations in terms of variables averaged across the thin gap. It is also assumed that the displaced viscoelastic fluid wets the wall, leaving films of uniform thickness behind. From this analysis, the elasticity of fluid and thus, the effects of the normal stress differences and shear stress thinning on the thin liquid film can be determined. In this paper the gas-assisted displacement of a non-Newtonian fluid described by the third-order fluid model confined in a Hele–Shaw cell was studied analytically using a singular perturbation method. The theoretical results obtained from this analysis are in agreement with experimental data [12] and theoretical results [11]. 2. Fluid model By using arguments of modern continuum thermodynamics, the model for a fluid of third-order has been derived by Bird et al. [20]. They show that the constitutive relation for the stress tensor for an incompressible homogeneous fluid of third-order has the form: h i s ¼
gcð1Þ þ a1 cð2Þ þ a2 fcð1Þ  cð2Þ g þ b1 cð3Þ þ b2 fcð1Þ  cð2Þ þ cð2Þ  cð1Þ g þ b3 ðcð1Þ : cð2Þ Þcð1Þ ;

ð2aÞ

where cð1Þ ¼ c_ ¼ rv þ ðrvÞþ ; D þ c
fðrvÞ  cðnÞ þ cðnÞ  ðrvÞg; cðnþ1Þ ¼ Dt ðnÞ

ð2bÞ n P 1.

ð2cÞ

The second- and third-order fluids are, respectively, obtained from Eq. (2a) by considering all of the terms in the first line of Eq. (2a) and all of the terms shown explicitly in the equation. In the above equation g is the viscosity a1, a2, b1, b2 and b3 are the material constants, that can be called them retarded-motion constants, D/ Dt denotes the material time derivative, v denotes the velocity field, and c(1) and c(2) are the first and second rate-of-strain tensors. The above model shall not be considered as an approximation to a simple fluid in the sense of a retardation, but rather an exact model in the sense described by Bird et al. [20]. The Clausius–Duhem inequality holds and the specific Helmholtz free energy is a minimum when the fluid is locally at rest, which leads to the following restrictions on the material coefficients. g P 0;

a1 P 0;

b1 ¼ b2 ¼ 0;

b3 P 0;




pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi 24gb3 6 a1 þ a2 6 24gb3 .

ð3Þ

Therefore, the model of (2) reduces to h i s ¼
gcð1Þ þ a1 cð2Þ þ a2 fcð1Þ  cð2Þ g þ b3 ðcð1Þ : cð2Þ Þcð1Þ .

ð4Þ

3. Physical problem and mathematical formulation Consider the motion of a long bubble moving steadily in a horizontal channel of rectangular cross-section (Hele–Shaw cell) with a half distance, d, between the parallel plates, filled with a viscoelastic fluid of viscosity, g, as shown schematically in Fig. 1. In this analysis, it is assumed that the displaced viscoelastic fluid leaves a thin film of uniform thickness b on the wall and the displacing fluid is invicid. The third-order fluid model is used to describe the viscoelastic fluid [20]. Total stress is given as follows: h i s ¼
gcð1Þ þ a1 cð2Þ þ a2 fcð1Þ  cð2Þ g þ b3 ðcð1Þ : cð2Þ Þcð1Þ ; ð5aÞ

1392

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

Fig. 1. Schematic diagram of gas-assisted displacement for different region (II: transition region, III: capillary static region).

where D þ c
fðrvÞ  cð1Þ þ cð1Þ  ðrvÞg; Dt ð1Þ Dcð1Þ ocð1Þ ocð1Þ ocð1Þ ¼ þu þv . Dt ot ox oy

cð2Þ ¼

ð5bÞ ð5cÞ

The gas-assisted displacement process can be considered as a singular perturbation problem. This approach for the case of a Newtonian fluid was first derived using a perturbation expansion in Ca1/3 [5]. Later the approach of Bretherton was expanded in both Ca1/3 and d = b/d and made some improvement in determining the film thickness and the pressure drop at the nose of the bubble [6,8]. Recently, the same method was applied to a viscoelastic fluid modeled by an Oldroyd-B constitutive equation [11]. In the perturbation method it is assumed that the bubble is moving with a sufficiently small value of the capillary number Ca and the interface is almost parallel to the wall of the channel. Therefore, the subsequent analysis employs the continuity equation and equation of motion for two-dimensional flow. Inner and outer expansions are developed in terms of small parameters d (d = De/Ca1/3), d2 and Ca1/3 for a non-Newtonian fluid. The method of matched asymptotic expansions is used to match the inner and outer solutions [21,22]. In this analysis, the origin of the frame of reference is taken to be the nose of the bubble moving in the positive x-direction. Consider the motion of a gas bubble into an incompressible Newtonian or non-Newtonian liquid as shown schematically in Fig. 1. For convenience, velocities are non-dimensionalized by the uniform velocity ub, the transverse and axial coordinates by the characteristic length d, and the pressure by r/d. The characteristic length d is taken to be half of the distance between the parallel plates. The equations of continuity and motion in the dimensionless form are given as follows: o
u o
v þ ¼ 0; ð6Þ o
x o
y    2   3    o
u o
u o
p o
o
u o3
u o3
u o3
u u o2
u ReCa
u þ
v þ þ De

u 2 3 þ ¼
þ Ca

v 2 2 þ 3 o
x o
y o
x o
x2 o
y 2 o
x o
xo
y 2 o
x o
y o
y   2    2   2  o
v o
o
u o
u o
u o
u u þ þ ð4 þ 16q1 Þ þ ð1 þ 4q1 Þ o
y o
y 2 o
x o
y o
x2 o
y o
x ( "    2  2 #   2  2 2 o
u o
u o
v o
u o
u o
u u 5=3 2 o
þ 16 24 þ4 þ8 þ q2 Ca De 2 o
x o
y o
y o
x o
y o
x o
y o
x " #   2   2  2  2    2 ) o
v o
v o2
o
u o
u o
v o
u o
v o
v u þ6 þ4 þ ; þ8 þ 2 4 o
y o
x o
y o
y o
y o
y o
y o
y o
x o
y 2 ð7Þ

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

 2     2    2  o
v o
v o
p o
v o2
v o3
u o3 v o
u o
u o
v o
u ReCa
u þ
v þ 2
De
u 2 þ v 3 þ
¼
þ Ca 2 o
x o
y o
x o
y o
y oy o
x o
xo
y o
x o
y 2 o
x o
y   2    2    2   2   o
o
v o
u o
v o
u o
v o
v o
u u
10
2
2
2 2 o
y o
y o
x o
x o
x o
x2 o
xo
y o
xo
y   2    2  o
v o
v o
u o
u þ ð2 þ 16q1 Þ þ 4q1 o
y o
y o
y 2 o
y 2 ( "   #     2 "  2 2 2 2 o
u o
v o2
u o
u o
u v 5=3 2 o
8 þ 24 þ6 þ q2 Ca De þ 4 2 o
x o
x o
x o
y o
xo
y o
y  2   #     2    2  o
v o
u o
v o
u o
u o
u o
v o
v þ 16 þ4 þ 8 þ8 2 o
y o
x o
y o
y o
x o
y o
x o
xo
y   2 )  2   o
v o
u o
v o
u . þ4 þ8 2 o
y o
y o
y 2 o
y

1393



ð8Þ

q1 is taken to be equal to
w2/w1 and w1 and w2 are the first and second normal stress coefficients, q2 is 4b3 gCa
2=3 =w21 and
p is the pressure.
u and
v are the dimensionless velocity components in
x and
y directions, respectively. The
y axis is taken normal to the channel plates with the origin at the mid-plane. Thus
y has the value of ±1 at the solid boundaries. The origin is assumed to be located at the nose of the bubble and consequently the velocity is independent of time with respect to this frame of reference. For obtaining Eqs. (6)–(8), the dimensionless variables are defined as follows:
x ¼ x=d;


y ¼ y=d;

H ¼ H =d;


u ¼ u=ub ;


v ¼ v=ub ;


p ¼ p=ðr=dÞ.

The above formulation is written with the assumption that the y derivative of the leading shear stress is the only significant force in the lubrication approximation. By the following analysis, the order of magnitude of the correction due to a normal stress gradient can be estimated. As can be seen from Eqs. (7) and (8) the forces due to inertia are negligible as Reynolds number goes to zero. Also the gravitational force has a negligible effect on the flow. Therefore, non-dimensionalization of the above equations and boundary conditions give rise to the following dimensionless parameters:   gub w1 ub Ca ¼ ; De ¼ K; K ¼ ; r 2g d where De denotes the Deborah number and w1 the first normal stress coefficient. The definition of w1 in terms of stresses will be given later. The Deborah number may be interpreted as the ratio of magnitude of the elastic forces to that of viscous forces, while K appears in the definition is a typical shear rate. Re = (dqub/g) was taken to be less than unity. In this analysis, the no-slip condition is applied on the solid boundaries. In order to establish the appropriate interfacial conditions, it is assumed that the viscosity of the gas within the bubble is negligible when compared with the viscosity of the fluid exterior to the bubble. It is also assumed that the displaced non-Newtonian liquid wets the walls of the channel and leaves a film of uniform thickness on the walls. In this analysis, the problem is solved using a perturbation method in expansions of the capillary and Deborah numbers in powers of Ca1/3, d (d = De/Ca1/3) and d2 which correspond to large surface tension and weak elasticity. The kinematic boundary condition at the interface is given by dy v ¼ . dx u

ð9Þ

The tangential stress boundary condition at the interface is independent of frame of reference and is given by ðp  nÞ  t ¼ 0. ð10Þ The difference in the normal stress on the two sides of the interface is balanced by surface tension and is expressed as

1394

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

ðp  nÞ  n ¼ p0
r=RP ;

ð11Þ

where p is the total stress tensor and is given by pij = pdij + sij, and n and t are unit vectors normal and tangent to the interface, respectively and subscript x and y denote their components. Also 1/RP is the principal radius of curvature of the interface and p0 is the pressure within the gas bubble. 3.1. Capillary static region For the capillary static region substituting expressions for p, t and n into Eqs. (10) and (11) and writing in dimensionless form gives: 
2 o
h oh o
h

sxy ¼ 0 at
y ¼

hðxÞ; þ
sxy

syy o
x o
x o
x " " 
2 #
1 " 
2 
 # 
2 #
3=2 oh oh oh o2
h oh
syy þ
sxx þ 2
sxy ¼
2 1 þ D
p
Ca 1 þ o
x o
x o
x o
x o
x

ð12Þ


sxx

at
y ¼

hðxÞ;

ð13Þ

where D
p ¼ p0
p,
u ¼
1;
v ¼ 0 at
y ¼ 1.

ð14Þ

In order to accomplish the matching between the capillary static region and the transition region these boundary conditions are augmented in later section. In order to do appropriate scaling, the flow region can be divided into four regions as shown in Fig. 1. The flow regions I, II, III and IV represent the thin film region, the meniscus front region, and the parallel or tubular flow region, respectively. The meniscus front region can be divided into two subregions: for small Ca, a capillary static region (III) where the shape of the interface is almost circular and the transition region (II) where the circular shape of the interface is deformed by viscous forces and shear thinning effects. The hydrostatic, nearly circular interface in the capillary static region is due to the dominance of pressure and interfacial tension in that region. The presence of the transition region allows matching of the solution in the capillary static region to the solution in the constant film thickness region. As will be seen later, the transition region is important since the interface is deformed by the shear thinning effects. In order to obtain the leading order terms in the outer expansion the Ca and d are taken to be zero (d ! 0, Ca ! 0) and Eqs. (7) and (8) become o
p00 ¼ 0 and o
x

o
p00 ¼ 0; o
y

ð15Þ

where the superscript (0) denotes the zeroth-order approximation in the outer expansion. Eq. (15) implies that the leading order approximation to the pressure
p00 ð
x;
y Þ is constant. In order to determine the value of this constant, the interfacial boundary conditions have to be examined when d ! 0 and Ca ! 0. As can be seen in Eq. (13), the viscous force term goes to zero when d ! 0 and Ca ! 0 and the gas bubble fills the entire channel. This implies that k ! 1 and the radius of curvature will be equal to the half distance between parallel plates. The lubrication approximation requires that the fluid move at very slow velocity and relative to the half distance between parallel plates, the liquid film thickness be small enough to assume that the interface is almost parallel to the wall of the channel and thus, the fluid motion in this region can be treated as a two-dimensional flow. When d ! 0 and Ca ! 0, Eq. (13) becomes 2 ! 3
3=2 00 00 2 o2
o
h h 4 00 5 D
p ¼
2 1þ . ð16Þ o
x o
x In order to obtain the static meniscus, the above equation is integrated subject to the conditions: 00
h
x ! 1

and

00
h ¼ 0 at
x ¼ 0.

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

1395

The solution of Eq. (16) is 1 00
h ð
xÞ ¼ 00 ½1
ðD
p00
x þ 1Þ2 1=2 . D
p

ð17Þ

Since the bubble fills the entire channel and the pressure is equal to the radius of curvature, the pressure must, therefore, be equal to
1. As can be seen in later section a higher-order expansion is required for the outer solution to match appropriate terms in the inner expansion. This can be done by substituting a higher-order expansion of the variables in powers of Ca1/3 and d into Eq. (13). The expansions of all unknown quantities in powers of d and Ca1/3 can be performed as follows: 1 X i;j pðx; yÞ ¼ pðx; yÞ di Cað1=3Þj ; i;j¼0

vðx; yÞ ¼

1 X

vðx; yÞi;j di Cað1=3Þj ;

i;j¼0

uðx; yÞ ¼

1 X

i;j

uðx; yÞ di Cað1=3Þj ;

i;j¼0

hðx; yÞ ¼

1 X

i;j

hðx; yÞ di Cað1=3Þj .

i;j¼0

It can be seen that the viscous terms will not make any contribution until powers of O(Ca) and d then the higher-order expansions of the variables in powers of Ca1/3 and d substituted into Eq. (13). These mean that D
p10 , D
p11 , D
p12 , D
p20 , D
p21 and D
p22 are constant and equal to perturbations in the radius of curvature. D
p10 , D
p11 , D
p12 , D
p20 , D
p21 and D
p22 are obtained for the expansions of O(dCa0), O(dCa1/3), O(dCa2/3), O(d2Ca0), O(d2Ca1/3) and O(d2Ca2/3) and they are given as follows: 10 2
1=2
p10
x½1
ð
x þ 1Þ  ; h ðxÞ ¼ D
11

ð18Þ

2
1=2

;

ð19Þ

2
1=2

;

ð20Þ

2
1=2

;

ð21Þ

21 2
1=2
h ðxÞ ¼ D
p21
x½1
ð
x þ 1Þ  ;

ð22Þ


h ðxÞ ¼ D
p11
x½1
ð
x þ 1Þ  12


h ðxÞ ¼ D
p
x½1
ð
x þ 1Þ  12

20


p20
x½1
ð
x þ 1Þ  h ðxÞ ¼ D
22

2
1=2


h ðxÞ ¼ D
p22
x½1
ð
x þ 1Þ 

ð23Þ

.

Taylor series expansions of Eqs. (18)–(23) about
x !
1 are given, respectively, as follows: 10 2 3
p10 þ D
p10 ð
x þ 1Þ
Dp10 ð
x þ 1Þ þ Oðð
x þ 1Þ Þ; h ð
xÞ ¼
D
11

2
h ð
xÞ ¼
D
p11 þ D
p11 ð
x þ 1Þ þ Oðð
x þ 1Þ Þ; 12

12

20

ð26Þ

2

ð27Þ

2

ð28Þ

2

ð29Þ


h ð
xÞ ¼
D
p20 þ D
p20 ð
x þ 1Þ þ Oðð
x þ 1Þ Þ; 21


h ð
xÞ ¼
D
p þ D
p ð
x þ 1Þ þ Oðð
x þ 1Þ Þ; 21

22

21

ð25Þ

2


h ð
xÞ ¼
D
p þ D
p ð
x þ 1Þ þ Oðð
x þ 1Þ Þ; 12

ð24Þ


h ð
xÞ ¼
D
p22 þ D
p22 ð
x þ 1Þ þ Oðð
x þ 1Þ Þ. 3.2. Transition region

In the transition region, the problem has to be rescaled in order to apply the lubrication approximation. In that region, the pressure, viscous force, interfacial tension and shear thinning behavior are all important. The inner variables may be rescaled in a manner similar to the Newtonian case as follows:

1396

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

~ u¼
u;

~v ¼


v ; Ca1=3

~p ¼
p;

~x ¼

ð
x þ lÞ ; Ca1=3


~h ¼ ð1
hÞ ; Ca2=3

~y ¼

ð
y þ 1Þ . Ca2=3

These rescaled variables reveal that the relationships between the inner and outer stresses are: ~sxx ¼ Ca1=3
sxx ;

~sxy ¼ Ca2=3
sxy ;

and ~syy ¼ Ca1=3
syy .

For the sake of simplicity from this point onwards, the tilde on the transition region variables shall be omitted and the letters x, y, u, v, p, h, sxx, sxy and syy be used to denote the dimensionless transition region variables. The origin of the coordinate system for the transition layer is moved to a new position ~x ¼ 0 as shown in Fig. 1. Equations of continuity and motion (Eqs. (7) and (8)) with rescaled variables are given by ou ov þ ¼ 0; ox oy op osxx osyx ¼
Ca2=3
; ox ox oy op osxy osyy ¼
Ca2=3
Ca2=3 . oy ox oy

ð30Þ ð31Þ ð32Þ

As can be seen from the rescaled momentum Eqs. (31), (32) as Ca ! 0, the leading viscous stress is the y derivative of shear stress that is the only significant force in the lubrication approximation. The stress tensor components in terms of the rescaled variables are given by ( "    2  2  2 #) 2 ou De o2 u o2 u ou 2 ou ou 2 ou sxx ¼
2 þ 1=3 2u 2 þ 2v
4
2=3
q1 8 þ 2=3 ox Ca ox oxoy ox oy ox oy Ca Ca  2  2  2 #  2  " De ou ou ou ov þ2 þ 4Ca2=3
2q2 4Ca2=3 ; ð33Þ 1=3 ox ox oy oy Ca  2       ou De ou o2 u ou ov ou ov þv 2
2 syx ¼
þ 1=3 u
2Ca2=3 oy Ca oxoy oy oy oy ox ox  2  2  2 #  2  " De ou ou ou ov þ2 þ 4Ca2=3
q2 4Ca2=3 ; ð34Þ oy ox oy oy Ca1=3 ( "     2  2  2 #) ov De o2 v o2 v ov ou ov 2 ou ov þ 2v 2
2 syy ¼
2 þ 1=3 2u þ q1 þ8
4 oy Ca oxoy oy ox oy oy oy Ca2=3 oy  2  2  2 #  2  " De ov ou 2=3 ou 2=3 ov þ2 þ 4Ca 4Ca
2q2 . ð35Þ oy ox oy oy Ca1=3 Note that stresses sxx and syy in Eqs. (33) and (35) contain elastic terms of a potentially large size of dCa
2/3 relative to the viscous stresses in the same directions. However, these ultimately do not dominate the overall stresses for small elasticity (d  1) since the normal stress terms are ultimately down-weighted by a capillary number in power of Ca2/3 relative to shear stress in the rescaled momentum equations (26), (27), and as expressed earlier the leading viscous stress is a shear stress sxy. The boundary conditions at the interface (y = h(x)) for the inner solution are given as follows: Total tangential force has to be equal to zero at the interface (y = h(x)) and with the inner variable it is given by    2   oh oh oh
2=3 syy þ Ca syx
sxx
syx ¼ 0 at y ¼ hðxÞ. ð36Þ ox ox ox The difference in the normal stress on the two sides of the interface is balanced by surface tension and for the inner solution it is expressed as follows:

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

" Dp
1 þ Ca

2=3



oh ox

2 #
1 " Ca

2=3

"  2 #
3=2 o2 h 2=3 oh ¼ 2 1 þ Ca ox ox

syy þ Ca

4=3

1397

 2   # oh 2=3 oh sxx
2Ca sxy ox ox

at y ¼ hðxÞ.

ð37Þ

The kinematic boundary condition at the interface is given by oh v ¼ ox u

at y ¼ hðxÞ.

ð38Þ

The boundary conditions for velocities: u ¼
1;

v ¼ 0 at y ¼ 0.

ð39Þ

Any additional conditions require besides the above conditions can be obtained from the matching conditions somewhere in the transition region. 3.3. The effects of elastic terms on the residual liquid film thickness To establish the magnitude of the leading order elastic term in each region, the stress equations are substituted into the equations of motion in the capillary static and transition regions. For small elasticity and large surface tension, the leading elastic terms in the capillary static and transition regions are determined to be of order CaDe and De/Ca1/3, respectively. The above analysis indicates that the addition of fluid elasticity has the grater impact on the shape of the free surface profile in the transition region for small capillary number. In the capillary static region, the weak elastic force, O(CaDe), has little effect on the hydrostatic circular interface shape in the presence of the strong surface tension. Therefore, the leading viscoelastic corrections to the pressure jump and thin film thickness in region-I are a result of the modified flow in the transition region. From the inner dependent and independent variables in the transition region, dimensional y- and x-axes must be scaled by dCa2/3 and dCa1/3, respectively and clearly b is of the same order as y in the inner region and thus, must be scaled by dCa2/3. d in terms of the length scales can be written as: d = De/Ca1/3 = (w1/2g)(ub/dCa1/3). Since dCa1/3 is the length scale for the axial direction in the transition region, the equation just given above indicates that d is the ratio between the characteristic time of the fluid (relaxation time), w1/2g, and time it takes for a fluid element to flow through the transition region (dCa1/3/ub). To the leading order, the residence time in the transition region is assumed to be much longer than the relaxation time of the fluid that corresponds to small Deborah number. Therefore, it can be said that the polymeric fluid shows only minor qualitative difference from a Newtonian fluid since the motions keep the polymer molecules more or less in their equilibrium configurations. In such a case, Newtonian fluid behavior is obtained in the limit as De ! 0. The effect of the finite relaxation time which is long enough to disturb the local flow and stress field away from their Newtonian values can be seen by examining the dynamics at O(d) and O(d2). The leading order shear stress in the transition region is obtained by neglecting time derivative in Eq. (5c). This quasi-steady assumption, justified because the entire coordinate system is advecting with the bubble leads to the determination of the leading order shear stress from the local shear rate. As stated earlier, at O(d) the shear stress of the fluid in transition region shows only minor qualitative difference from that of a Newtonian fluid due to the time it takes for a fluid element to be advected to a zone of a different shear rate being comparable to the relaxation time of the fluid. 3.4. Obtaining the governing equations As Ca ! 0 and d ! 0, the leading order (zeroth-order) terms are obtained by substituting the expressions for stresses (Eqs. (33)–(35)) into Eqs. (30)–(32) and some of the boundary conditions namely the tangential stress, normal stress and kinematic boundary conditions.

1398

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

The equations of continuity and motion for the lowest order approximation are obtained as ou00 ov00 þ ¼ 0; ox oy

op00 o2 u00 ¼ ; ox oy 2

op00 ¼ 0. oy

ð40Þ

For zeroth-order approximation the tangential stress, normal stress and kinematic boundary condition are, respectively, given by following equation: ou00 o2 h00 oh00 ¼ 0; Dp ¼ ¼ v00 ; u00 2 oy ox ox

at y ¼ h00 ðxÞ.

ð41Þ

These equations are given elsewhere for a Newtonian fluid [6,5] and for a power-law fluid [23]. The following well-known third-order differential equation is obtained from Eqs. (40) and (41) and the related boundary conditions. That is used to determine the interface profile and thus, the fraction of liquid deposited on the wall of the geometry: o3 h00 3ðh00
b00 Þ ¼ ; 3 ox3 ðh00 Þ

ð42Þ

where b00 is residual liquid film thickness at the constant film region in gas-assisted displacement process. The additional boundary and matching conditions to which Eq. (42) is subject will be discussed in more detail later. Similarly as Ca ! 0 and d1, the leading order O(d1Ca0) terms are obtained by substituting the expressions for stresses (Eqs. (33)–(35)) into Eqs. (30)–(32) and using the related conditions. For the expansion of O(d1Ca0) the equations of continuity and motion are given by ou10 ov10 þ ¼ 0; ox oy

ð43Þ

   00  2 00  2 00 3 00 o2 u10 op10 ou00 o2 u00 ou ou 00 o u 00 o u þ u þ v ¼ þ Þ
ð1
4q ; 1 2 3 2 oy ox oyox oy ox oy oy oyox op10 ¼ 0. oy

ð44Þ ð45Þ

The above equations are obtained by considering variables with expansion of O(dCa0) that corresponds to a non-Newtonian fluid described by the second-order fluid model. The three-interface boundary conditions namely the tangential stress, normal stress and kinematic boundary conditions for the expansion of O(d1Ca0) are, respectively, given by    00 2 00 2 00 2 00 ou10 o2 u00 ou00 ou00 ou oh 00 o u 00 o u ¼
h10 þ v þ u þ 2 þ 2 2 2 oy oy oyox oy oy ox oy ox o2 h10 at y ¼ h00 ðxÞ; ox2     Z h00 oh10 o 1 op00 o 2 ¼
u10 dy. ðh00 Þ h10 þ ox 2 ox ox 0 ox

Dp10 ¼

at y ¼ h00 ðxÞ

ð46Þ ð47Þ ð48Þ

The boundary conditions for velocities are u10 ¼ 0 and v10 ¼ 0 at y ¼ 0.

ð49Þ 10

As in a power-law fluid case, a third-order differential equation in h may be obtained from Eqs. (43), (44), (47), (48), and the associated boundary conditions (Eqs. (46) and (49)) [13]:

00 2 b o3 h10 6h10 9b00 h10 3b10 h00 b00 00 X ¼
3 þ 4
3 þ ð3
0:9q1 Þ 3
ð18
9q1 Þ 4 hX þ ð18
8:1q1 Þ 5 h00 X . 3 00 00 00 00 00 ox h00 h h h h h ð50Þ

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

1399

The integration constant is obtained from the boundary condition subject to h10 = b10 as x !
1. The above third-order differential equation governs the gas-assisted displacement of a viscoelastic fluid described by the second-order fluid model. The residual liquid film thickness of the second-order fluid will be obtained by the numerical solution of Eq. (50). Now the third-order differential equation valid for a viscoelastic fluid described by the third-order fluid model can be obtained by considering the variables with expansions of O(d2Ca0), O(d2Ca1/3) and O(d2Ca2/3). Similar to the above derivation, as Ca ! 0, and d2, the leading order O(d2Ca0) terms are obtained by substituting the expressions for stresses (Eqs. (33)–(35)) into Eqs. (30)–(32) and using the related conditions. For the expansion of O(d2Ca0) the equations of continuity and motion are given by ou20 ov20 þ ¼0 ox oy

ð51Þ

    3 10 3 10 3 10 3 00 o2 u20 op20 ov00 o2 u10 ov10 o2 u00 00 o u 10 o u 00 o u 10 o u þ u þ u þ v ¼ þ v

oy 2 ox oy 2 ox oy 2 ox oy 3 oy 3 oy oy 2 oy oy 2  2 00  00 2  00  2 10   10  2 00  ou ou ou ou ou ou
ð1
4q1 Þ ; þ
6q2 oy oyox oy oyox oy 2 oy op20 ¼ 0. oy For the expansion of O(d2Ca0) the boundary conditions at the interface are given by   2 00 2 10 2 10 2 00 2 10 2 00 ou20 ou00 ou10 20 o u 10 o u 00 o u 10 o u 00 o u 10 o u ¼
h þu þv
h þu þv þ2 oy oy 2 oy 2 oyox oyox oy 2 oy 2 oy ox         2 2 00 10 ou10 ou00 ou10 ou00 oh ou00 oh ou00 ou00 þ2
2q2 þ2 at y ¼ h00 ðxÞ; þ4 oy ox oy oy ox oy ox oy oy o2 h20 at y ¼ h00 ðxÞ; ox2     Z h00  20  oh20 o 1 op00 00 2 20 ou10 10 o2 u00 10 2 ou ¼
h þ h h þ þ h dy. ox 2 ox ox ox oxoy ox 0

Dp20 ¼

ð52Þ ð53Þ

ð54Þ ð55Þ ð56Þ

The boundary conditions for velocities are u20 ¼ 0 and v20 ¼ 0

at y ¼ 0.

ð57Þ

The third-order differential equation in h20 may be obtained from Eqs. (51), (52), (55), (56) and the related boundary conditions (Eqs. (54) and (57)). Eqs. (42), (50) and (58) can be expressed in translation form since they are autonomous equations. The following transformation is introduced: X ¼

x þ x0 b00

and

H i;j ðxÞ ¼

hi;j ðxÞ b00

Eqs. (59)–(61) are obtained by introducing the above new variables into Eqs. (42), (50) and (58). Since Eqs. (59)–(61) are third-order equations; the solution of each of these equations will contain three constants. The solution of Eq. (59) as X !
1 is given by H 00 ðX Þ ¼ 1 þ expð31=3 X Þ
ð3=7Þ expð2x31=3 X Þ þ    . It is clearly seen from the above equation H00 ! 1 as X !
1 (for more details see [5]). A solution in the transition layer has to be matched to both a solution in thin film and that in the capillary static region. Matching of the inner solution to that for the thin film is accomplished by imposing the boundary condition of H00 ! 1 as X !
1. The method of matched asymptotic expansions requires that the inner solution as X ! 1 match to the outer solution as
x !
‘. Therefore, it is necessary to find asymptotic expansion of Eq (59) as X ! 1. In

1400

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

the limit as Ca ! 0, the bubble fills the entire region of the geometry. The angle of contact at the wall must be zero; thus RP = ± 1 and the capillary static profile is circular. Thereby, it can be expected that Hi,j(X) eventually becomes the following quadratic form for large X. In other words, it is possible to show that as X approaches infinity, H00, H10 and H20 have the following quadratic forms (Eqs. (62)–(64)).

2

2

2 b00 h10 o3 h20 6h20 9b00 h20 3b20 9 h10 9 b10 h10 ¼
3 þ 4
3 þ 4
9 5 þ

 2 h00 2 h00 4 ox3 h00 h00 h00 h00 "   10   1 231 63 h 135 1065 b00 h10 00 q1 þ q
3hx ð1 þ q1 Þ 3
þ þ

00 5

00 4 20 4 10 1 40 h h00 h     00 2 10   10 h 9 27 b 15 15 b 63 1 q
q1 þ q1 þ 45 4


00 6 þ

00 4
00 20 1 8 2 8 10 h h h   00   00 2   00 3 b b 144 2 189 b 144 2 189 63 q þ q þ 171 5
q1 þ q1 þ 171 6 þ q1 þ 45 7 þ 00 00 5 1 2 1 5 2 10 h h00 h #   00   2   36 2 231 b 36 2 117 27 b00 99 63 b00 b10 q þ q þ 15 4 þ q þ q þ q þ


00 5


00 5 5 1 10 1 5 1 5 1 2 10 1 8 h h h00 "   

2 1 219564 2 10461 1491 1 304 2 q þ q
q1 þ 12 h00
3h00

00 3 þ 60q1 þ xx x 1575 1 3780 1 20 5 h00 h     00

00 2 b00 170142 2 377685 10497 b 608 2 q1
q1
q

120q þ þ 24 hx  2

00 4 1 350 2100 40 5 1 h h00   2   3 # 98214 2 80253 11808 b00 75078 2 12835 4293 b00 q
q
q
q
þ

00 5


00 6 175 1 420 1 40 350 1 2100 1 40 h h    00   

2 304 152 b q þ 4 þ 30q21 þ q þ 6 00
3 h00
20q21 þ xx 15 1 5 1 h "    00

00 2 1922895 2 93465 1827 1 7252869 2 404499 28587 b q1 þ q þ q1
q1


3 hx

00 5

00 4 þ 700 560 1 5 700 560 20 h h   00 2   00 3 40084785 2 1716489 87675 b 18088515 2 413451 27477 b q1 þ q þ q1
q1
þ


00 7

00 6 þ 3500 1680 1 40 3500 1680 40 h h #

00 2

 2 b 54 15 h10 45 b00 h10 15 b10 135 b00 h10 45 b00 b10
q2  7 þ
þ þ





00 7




 5 2 h00 5 2 h00 6 4 h00 5 8 8 h00 6 h00 h "    00   00 2 #

10  b 3 1 9 b 27 20
q1 3
42 þ q1 4 þ 18
q1 5
3 hx 00 00 10 5 10 h00 h h "     00   00 2

00 3 b 15 9 1 225 207 b 225 þ q1  5 þ þ q1  6
þ 27q1 7
3 hx
00 00 4 8 8 16 4 h h00 h

00 2   00 3 # b b 135 243 162 1 486 b00 486 þ q1  8 þ q2  4
q2 5 þ q2  6 ; þ ð58Þ 00 00 00 4 16 5 5 5 h h h h00 H 00 XXX ¼

3ðH 00
1Þ ðH 00 Þ

3

;

ð59Þ

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

H 10 XXX ¼


6H 10

9H 10

3ðb10 =b00 Þ

ðH 00 Þ

ðH 00 Þ

ðH 00 Þ

þ 3

þ ð18
8:1qÞ H 20 XXX


4

H 00 X b ðH 00 Þ5 00

;

3

þ ð3
0:9qÞ

H 00 X b00 ðH 00 Þ


ð18
9qÞ 3

1401

H 00 X b00 ðH 00 Þ

4

ð60Þ

2

10 2 H 6H 20 9H 20 3b20 =b00 9 H 10 9 b10 =b00 H 10 ¼


þ
þ
9 þ 3 00 4

00 3



00 4

00 5 2 H 00 4 2 H 00 H H H h "     10 27 1 693 189 H 405 3195 H 10 1 þ q q q
H 00 ð Þ
þ þ þ

00 5

00 3

00 4 1 1 1 X 00 00 00 2 20 4 b H 10 40 b H b H       10 00 10 27 81 H 45 45 b =b 189 1 q
q þ q þ 135 2




 þ


4 00 20 1 8 b00 H 00 6 2 1 8 b00 H 00 4 10 1 b H 00     432 2 567 1 432 2 567 1 q1 þ q1 þ 513 2

q q þ
þ þ 513 



00 6 1 1 2 5 00 00 5 2 5 2 b b H 00 H     189 1 108 2 693 1 q þ 135 2

q þ q þ 45 þ

4 7
10 1 5 1 10 1 b00 H 00 b00 H 00 #     108 2 351 81 1 297 189 b10 =b00 q1 þ q1 þ q þ
þ



 5 5 2 b00 H 00 5 10 1 8 b00 H 00 5 "   

00 2 1 219564 2 10461 4473 1 912 2 q q q
H 00 þ
þ 180q þ þ 36 HX

00 2 00 3 1 1 XX 1 1 525 1260 20 5 H 00 b H    

2 1 510426 2 377685 31491 1 1824 2 q
q
q1 þ 24 H 00


00 2

00 2 00 4
360q1 þ X 350 1 700 1 40 5 H b H     294642 2 80253 35424 1 22521 2 12835 12879 q
q
q
q
þ

00 2 00 5
175 1 140 1 40 350 1 700 1 40 b H #      

00 2 1 304 456 1 2 2 q q  2


60q þ þ 12 þ 90q þ þ 18
H 6 XX 1 1 5 1 5 1 H 00 b00 H 00 " 

2 5768685 2 280395 5481 1 q q
þ þ
H 00

00 2 00 4 X 1 700 560 1 5 b H     21758607 2 1213497 85761 1 24050871 2 1716489 263025 q1
q1
q1 þ q1 þ þ

00 2 00 5 þ
700 560 20 700 560 40 b H   1 54265545 2 413451 82431 1 162 1 q1
q1
q2 3

 2

 þ

00 2 00 7 þ 7 00 00 00 6 3500 560 40 5 b b b H H H 00 # 45 H 10 135 H 10 45 b10 =b00 405 H 10 135 b10 =b00 þ








 þ



 þ








 2 b00 2 H 00 5 2 b00 2 H 00 6 4 b00 2 H 00 5 8 b00 2 H 00 7 8 b00 2 H 00 6 " #     

10  9 1 27 1 81 1 60
q1
HX

3
126 þ q1 00 00 4 þ 54
q1 00 00 5 10 5 10 b H b00 H 00 b H "    

3 45 27 1 675 621 1 þ q1 3

þ q 1 3


H 00
5 þ 6 X 00 00 00 4 8 8 16 b b H H 00 #     675 1 405 72 1 162 1 þ 81q1 3 7 þ þ q1  3

q2  2


þ  4 8 00 00 00 00 4 4 16 5 b b b h H 00 H 00 486 1 486 1 q 

q 


ð61Þ  þ  . 5 2 b00 2 H 00 5 5 2 b00 2 H 00 6

1402

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

The constants of the quadratic equations are determined from numerical integration of Eqs. (59)–(61). As X ! 1, the solutions of Eqs. (59)–(61) can be, respectively, expressed as follows: H 00 ¼ ðAN =2ÞX 2 þ BN X þ C N ;

ð62Þ

H 10 ¼ ðAS =2ÞX 2 þ BS X þ C S ;

ð63Þ

H 20 ¼ ðAT =2ÞX 2 þ BT X þ C T ;

ð64Þ

where N indicates Newtonian (the first-order fluid model), S the second-order and T the third-order fluid models. The inner solutions and outer solutions have now been determined separately. The inner and outer solutions have to match somewhere in the transition region. The asymptotic matching is accomplished at X ! 1 for the inner solution and at
x !
‘ for the outer solution. Since the relationship between the inner and outer dependent variables for the interface profile is given by h ¼ ð1

hÞCa
2=3 , Eq. (62) can be expressed in terms of the original variables as follows:       AN AN 2 00 2 00 1=3 AN 2=3



Ca x þ B x þ B x þ b C ð x þ ‘ Þ ð x þ ‘ Þ
Ca . ð65Þ h ¼1
0 N N 0 N 0 2b00 b00 2b00 00

The value of ‘ can be determined from the matching condition
h ! 1 as
x !
‘. This implies ‘ = 1 and the other equations for the outer solution may be expressed in a Taylor series expansion and matched with Eq. (65) to determine some of constants. However, this matching already exists in the literature [6] will not be repeated here. From the matching   b00 BN AN 2 00
x0 ¼
; b00 ¼ AN ; Dp01 ¼ 0; and Dp02 ¼ x þ B x þ b C ð66Þ N 0 N . 0 AN 2b00
10 is also accomplished at X ! 1 for the inner solution and at
x !
‘ for the The asymptotic matching for h outer solution. The inner solution, Eq. (63) is also expressed in terms of the original variables as follows:       AS AS 2 10 2 00 1=3 AS 2=3
x0 þ BS ð
x þ ‘Þ
Ca x0 þ BS x0 þ b C S . h ¼
ð
x þ ‘Þ
Ca ð67Þ 2b00 b00 2b00 If Eq. (24) is compared with Eq. (67) to powers of O(dCa0), it can be found that D
p10 must be equal to zero. Therefore, AS must also be equal to zero: D
p10 ¼ 0;

10
h ¼ 0;

and

AS ¼ 0.

Comparing Eq. (25) with Eq. (67) in terms of powers of O(dCa1/3) reveals that D
p11 ¼ 0

and

11
h ¼ 0.

From Eqs. (26) and (67) the following relationship is obtained up to powers of O(dCa2/3) as x ! 1 and
x !
1:   AS 2 00 D
p12 ¼ x þ B x þ b C ð68Þ S 0 S . 0 2b00 The constants of Eq. (63) can be determined by the numerical integration of Eq. (60). However, Eq. (60) has many solutions depending on the values of parameters b10/b00 and q1. As stated previously, the first constant of Eq. (63) has to be equal to zero; therefore, the solution has to be found which makes AS equal to zero. In fact the values of constants for Eq. (63) also depend on an interval of the integration used to integrate Eq. (60) numerically. However, the value of b10/b00 and q1 make AS zero does not depend on the interval of the integration. Thereby, the values of b10/b00 and q1 were determined from this limiting condition. 20 Similarly the asymptotic matching for
h is accomplished at X ! 1 for the inner solution and at
x !
‘ for the outer solution. The inner solution, Eq. (64), is also expressed in terms of the original variables as follows:

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408 20
h ¼


     AT AT 2 2 00 1=3 AT 2=3 x0 þ BT ð
x þ ‘Þ
Ca x 0 þ BT x 0 þ b C T . ð
x þ ‘Þ
Ca 2b00 b00 2b00

1403



ð69Þ

Comparison of Eq. (27) with Eq. (69) to powers of O(d2Ca0) makes D
p20 be zero. Therefore, AT must also be equal to zero: D
p20 ¼ 0;

20
h ¼ 0;

and

AT ¼ 0.

Comparing Eq. (28) with Eq. (69) in terms of powers of O(d2Ca1/3) reveals that
21 ¼ 0. D
p21 ¼ 0 and h From Eqs. (29) and (69) the following relationship is obtained up to powers of O(d2Ca2/3) as x ! 1 and
x !
1:   AT 2 00 22 D
p ¼ x0 þ BT x0 þ b C T . ð70Þ 2b00 In order to determine the value of p
22 the constants of Eq. (64) have to be determined. The constants of Eq. (64) can be determined by the numerical integration of Eq. (61). After determining the value of b10 that makes AS be zero from the numerical integration of Eqs. (59) and (60) simultaneously, the determined value of b10 is used in the numerical integration of Eq. (61). The differential equation (61) for H20 was simultaneously integrated with Eqs. (59) and (60) since some variables in Eq. (61) have to be determined from Eqs. (59) and (60). The differential equation (61) for H20 has also many solutions depending on the values of parameters b20/b00 and q2 as in Eq. (60). By applying matching conditions to the O(d2Ca0) capillary static region free surface profile equation (27) D
p20 ¼ 0;

20
h ¼ 0;

and

AT ¼ 0.

The value of b20/b00 was determined from the above condition. The sign of b20 indicates how the elasticity of a fluid affects the thin film thickness hydrodynamics of a viscoelastic fluid modeled by the third-order fluid approximation. The film thickness tends to increase as b20 has a plus sign, otherwise it tends to decrease. 4. Results and discussion Although the main purpose was to examine the effect of the fluid elasticity on the thin film hydrodynamics of a viscoelastic fluid described by the third-order fluid model when displaced by a long gas bubble in a Hele– Shaw cell,—since the obtained third-order governing differential equation from the use of the third-order fluid model contains some variables to be determined from the first and second-order fluid models—the first-, second- and third-order fluid approximations are examined together. The third-order differential equations (59)– (61) governing the gas-assisted displacement of the first-order (Newtonian), second-order and third-order fluids are obtained analytically. To obtain the viscoelastic free surface corrections, they must be integrated numerically and simultaneously. As stated previously, the first correction value of parameter b10/b00 that makes AS zero is determined from the numerical integration of Eqs. (59) and (60) simultaneously using a forth-order Runge–Kutta method by varying values of b10/b00 at a constant value of q1. The value of q1 has to be less than 0.5 and should be around 0.25 [20]. The values of AS as a function of the values of parameter b10/b00 at a constant value of q1 are shown graphically in Fig. 2. As can be seen in the figure the value of b10/b00 that makes AS zero is
0.1741. This result indicates that the film thickness of region-I has a small correction because of elasticity of a viscoelastic fluid that comes from the second-order fluid approximation. According to this result it may be said that the viscoelasticity of the fluid has a small effect on the fraction of liquid deposited on the wall at low Deborah numbers (De < 1). However, strictly speaking, our prediction states that the film thickness of a viscoelastic fluid tends to decrease as the fluid becomes more viscoelastic since the sign of b10 (the first viscoelastic correction term) is negative. On the other hand, the variation of AT values as a function of parameter b20 at a constant value of q2 is shown graphically in Fig. 3. The value of b20 where AT is zero for a constant value of q2 = 0.06 is clearly seen to be 0.80 in Fig. 3 and indicates that the second viscoelastic term affects the thin viscoelastic film

1404

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408 3

q 1 = 0.25

2.5

AS

2 1.5 1 0.5 0 -0.5 -0.6

-0.5

-0.4

-0.3

-0.2

-0.1

0

b10 Fig. 2. The values of the first constant of Eq. (63) AS obtained from the numerical integration of Eq. (60) for various values of b10 at the constant value of q1.

0.8 0.7

q 2 = 0.060

0.6

AT

0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 0

0.2

0.4

0.6

0.8

1

b20 Fig. 3. The values of the first constant of Eq. (64) AT obtained from the numerical integration of Eq. (61) for various values of b20 at the constant value of q2.

thickness in a positive way. By considering the total effect of the viscoelastic correction terms, the thin liquid film thickness of a viscoelastic fluid tends to decrease as the fluid becomes more viscoelastic since the numerator of the second correction term is much smaller that that of the first correction term in Eq. (72). Similar results were obtained experimentally by Huzyak and Koelling [12] and theoretically by Ro and Homsy [11]. Huzyak and Koelling experimentally examined the effect of fluid elasticity on the liquid film thickness by varying Deborah number from 0 to 5. They pointed out that at low Deborah numbers, De < 1, the value of reduced fractional coverage of viscoelastic fluids, mR = m/mN (m: fractional coverage of viscoelastic fluid, mN: fractional coverage of Newtonian fluid) is between 0.95 and 1.05, which means the film thickness of viscoelastic fluid is either slightly lower or slightly higher than that of Newtonian fluid depending on the reduced fractional coverage values taken to be either 0.95 or 1.05. The reduced fractional coverage, mR, less than unity indicates that the film thickness of a viscoelastic fluid is lower than that of a Newtonian fluid at an equivalent capillary number. Furthermore, the theoretical analysis of the gas-assisted displacement of a viscoelastic fluid modeled by an Oldroyd-B constitutive equation also showed that the liquid film thickness tends to decrease as the fluid becomes more viscoelastic [11]. In the present analysis, while the first viscoelastic correction term obtained from the second-order fluid approximation tends to decrease the liquid film thickness, the second viscoelatic correction term obtained from the third-order fluid approximation tends to increase the liquid film thickness as the fluid becomes more viscoelastic since the signs of b10 and b20 are negative and positive, respectively. The first and second viscoelastic correction terms come from the effect of normal stress and shear stress thinning, respectively. Although the

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

1405

0.1 b10 = -0.233

0.05

AS

0 -0.05 -0.1 -0.15 -1.5

-1

-0.5

0

0.5

1

1.5

q1

Fig. 4. The values of the first constant of Eq. (63) AS obtained from the numerical integration of Eq. (60) for various values of q1 at the constant value of b10.

second-order fluid model is able to predict the normal stress differences, it does not take into account the shear thinning and thickening phenomena that non-Newtonian fluids show. The third-order fluid model represents a further, although inconclusive, attempt toward a comprehensive description of the properties of viscoelastic fluids. In Fig. 4 the values of AS are obtained by varying values of q1 for the constant values of b10 from the numerical integration of Eq. (60). As can be seen in Fig. 4 the values of AS vary in a small interval (
0.045 to 0.072) while the values of q1 are varying from
1.0 to 1.0. As given previously, q1 equals
w2/w1 and the functions w1 and w2 are known as the first and second normal stress coefficients, respectively. The stress tensors depend only on the flow field; it can be said that the steady-state stresses in steady simple shear flow are functions only of the shear rate. The normal stress coefficients w1 and w2 are defined as follows: sxx
syy ¼
w1 ðcÞc2xy ;

syy
szz ¼
w2 ðcÞc2xy .

q1 in terms of the stress components is equal to
(syy
szz)/(sxx
syy) and szz = 0. As can be seen in Fig. 3 the value of q1 = 0.25 makes AS zero at the constant value of b10 =
0.23275. This indicates that stress in the axial direction is higher than that in transverse direction as expected. On the other hand, q2 is a material constant given by 4b3 gCa
2=3 =w21 . This constant was taken to be 0.06 that makes AT be zero at the value of b20 ffi 0.80 (Fig. 5). The value of q2 should be less than unity [20] and its sign has to be positive [24]. In Fig. 5 the values of AT as a function of q2 are depicted from the numerical integration of Eq. (61) at the constant value of b20. The following equation for the pressure drop at the static meniscus region can be written by considering the expansions of pressure in simple powers of d and Ca1/3 and matching conditions. The contribution terms from these expansions are

0.4 0.3

b 20 = 0.80

0.2

AT

0.1 0

-0.1 -0.2 -0.3 -0.4 0

0.02

0.04

0.06

0.08

0.1

0.12

q2 Fig. 5. The values of the first constant of Eq. (64) AT obtained from the numerical integration of Eq. (61) for various values of q2 at the constant value of b20.

1406

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

D
p ¼ D
p00 þ Ca2=3 D
p02 þ dCa2=3 D
p12 þ d2 Ca2=3 D
p22 ; where D
p00 ¼ 1=Rp, D
p02 , D
p12 and D
p22 are given by Eqs. (66), (68) and (70), respectively. From the above relationships the pressure at the nose of bubble is given by the numerical procedure as

r Dp ¼ 1 þ 3:8Ca2=3 þ 0:706DeCa1=3
0:805De2 . ð71Þ d Eq. (71) indicates that the first viscoelastic correction term of pressure comes from the second-order fluid approximation affects the pressure drop at the tip in a positive way while the second one comes from the third-order fluid approximation affects it in a negative way. In other words, while an increase in the normal stress will increase the pressure drop at the tip, an increase in the shear thinning will decrease it. By considering the total effect of the first and second viscoelastic correction terms of the pressure drop, the elasticity of a fluid tends to increase the pressure drop at the tip since the value of the parameter in front of the second correction term (De2) is much smaller than that in front of the first one (DeCa1/3). Thus, the pressure drop at the interface for a viscoelastic fluid displacement in a Hele–Shaw cell is slightly higher than that for Newtonian fluid displacement. It is also observed that the pressure drop slightly increases by increasing the value of q1 in the numerical integration of Eq. (60). As can be seen in Fig. 2 the value of b10/b00 that makes AS zero is equal to
0.1741. This reveals that the film thickness of region-I has a small correction term for the expansions of O(dCa0) and O(dCa1/3) even though the profile of the transition region does have a correction at these orders. The viscoelastic correction to the film thickness can be determined by selecting the O(d) and O(d2) free surface profile of transition region which matches the profiles of the constant film thickness region and capillary static region. Thus, the film thickness deposited on the walls can be expressed as follows: b ¼ 1:337Ca2=3
0:2327DeCa1=3 þ 0:825De2 þ Oðd3 Ca2=3 Þ. d

ð72Þ

Eq. (72) indicates that the effect of elasticity on the fraction of viscoelasticity liquid deposited on the wall during the gas-assisted displacement process is not considerable level since while in the right hand side of Eq. (72) the second term tends to decrease the film thickness, the third one tends to increase that of a viscoelastic fluid. From the comparison of the values of the first and second viscoelastic correction terms, it will be seen that the value of the second term (b20) obtained from the third-order fluid approximation is higher than that of the first viscoelastic correction term (b10) obtained from the second-order fluid approximation. This result indicates that the shear thinning is more effective than that of normal stress on the viscoelastic film thickness. As stated previously, while the first viscoelastic correction term tends to decrease the film thickness of viscoelastic fluid, the second one tends to increase the film thickness of viscoelastic fluid as the fluid become more viscoelastic. Although the value of the second viscoelastic correction term is higher than that of the first one, the film thickness of a viscoelastic fluid will be less than that of a Newtonian fluid at low capillary and Deborah numbers since the value of second term in right hand side of Eq. (72) is larger than that of the third one. The total effect of viscoelastic correction terms to the static profile in the capillary static region is much smaller than orders of magnitude considered in this paper. As mentioned earlier, the lubrication approximation is valid only for low capillary and Deborah numbers. The previous workers [15,16] reported that the instability of the interface for two-phase immiscible displacement of a Newtonian fluid in a planar geometry did not occur at low capillary number. One major assumption of present theory is that the displaced viscoelastic fluid wets the wall, leaving films of uniform thickness. This assumption rejects the instability of interface for two-phase immiscible displacement of a viscoelastic fluid confined in a Hele–Shaw cell. In present analysis, the residual liquid film thickness tends to decrease as the fluid becomes more viscoelastic at low capillary number (De < 1). The identical trend is observed in both experimental [12] and theoretical [11] studies. Even though both the theoretical analysis and the experimental data, in general, show that fraction of the viscoelastic liquid deposited on the walls, is less than that of Newtonian fluid at low Deborah numbers, a direct comparison of the experimental data with the theoretical results should be made with caution. Note that, when the origin of the reference frame is taken to be the nose of the bubble, the gas-assisted displacement problem in terms of determining the residual liquid film thickness on the wall becomes very

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

1407

similar to the free coating of a vertically withdrawn plate from the liquid reservoir. Therefore, the results of this theoretical analysis can be extended to the free coating process with small modifications. It must be kept in mind that the second-order fluid approximation cannot describe the relation between the non-Newtonian viscosity and strain faithfully. Thus, third-order fluid approximation is used to analyze the effect of fluid elasticity on the film thickness of a viscoelastic fluid deposited on the walls of Hele–Shaw cell. As a result, the present analysis qualitatively agrees with experimental data [12] and theoretical analysis [11] at low Deborah number. 5. Conclusions The second- and third-order fluid models are used to obtain the viscoelastic correction terms that describe ways in which Newtonian solutions are perturbed by weak viscoelasticity in the constant film region. It is concluded that the effect of the fluid elasticity on the film thickness of a viscoelastic confined in a planar geometry is relatively insignificant levels since the first viscoelastic correction term obtained from the second-order fluid approximation tends to decrease the film thickness of the viscoelastic fluid, the second correction term obtained from the third-order fluid approximation tends to increase the film thickness as the fluid becomes more viscoelastic. The value of the second viscoelastic correction term is higher than that of the first viscoelastic correction term. Therefore, the effect of shear thinning on the viscoelastic film thickness is more pronounced than that of the normal stress. An increase in shear stress will increase the film thickness of a fluid deposited on the wall. On the other hand, the normal stress force represents the restoring force of polymer, which acts to resist stretching in the flow direction; therefore, an increase in the normal stress force will decrease the liquid film thickness on the wall. By considering the total effect of fluid elasticity, it can be said that the effect of viscoelasticity is relatively insignificant level and a viscoelastic fluid behaves nearly identical to a Newtonian fluid at the low capillary and Deborah numbers. As mentioned earlier, at low capillary number (De < 1) the experimental study reported that a viscoelastic fluid exhibits a liquid film thickness identical to that of a Newtonian fluid at an equivalent capillary number. The first and second viscoelastic correction terms for the film thickness in the constant film region and the pressure drop in the capillary static region are obtained from the modified flow in the transition region. Since the total effect of viscoelastic correction on the static profile in the capillary static region is much smaller than orders of magnitude considered in this paper, the viscoelastic correction to the pressure jump condition has been obtained from the static solution of capillary static region. The theory presented in this paper reveals that the pressure drop at the tip decreases with increasing shear thinning of a viscoelastic fluid while it increases with increasing the normal stress force. One major assumption here, the displaced viscoelastic fluid wets the wall, leaving films of uniform thickness. This assumption rejects the instability of interface for two-phase immiscible displacement of a viscoelastic fluid confined in a Hele–Shaw cell. Thus, the presented theoretical analysis is valid only at low capillary and Deborah numbers. At low Deborah numbers the effects of elasticity of a viscoelastic fluid on the pressure drop at the interface are inconsiderable level also. This theoretical analysis in terms of determining the film thickness of viscoelastic liquid and pressure drop at the nose of bubble is in agreement with the theoretical analyses and experimental studies of some previous investigators. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

F. Fairbrother, A.E. Stubbs, J. Chem. Soc. 1 (1935) 527–529. G.I. Taylor, J. Fluid Mech. 10 (1961) 161–165. B.G. Cox, J. Fluid Mech. 14 (1962) 81–96. B.G. Cox, J. Fluid Mech. 20 (1964) 193–200. F.P. Bretherton, J. Fluid Mech. 10 (1961) 166–188. C.W. Park, G.M. Homsy, J. Fluid Mech. 139 (1984) 291–308. L.W. Schwart, H.W. Princen, A.D. Kiss, J. Fluid Mech. 172 (1986) 259–275. D.A. Reinelt, J. Fluid Mech. 183 (1987) 219–234. J. Ratulowskiand, H.C. Chang, Phys. Fluids A 1 (1989) 1642–1655. W.B. Kolb, R.L. Cerro, Chem. Eng. Sci. 46 (1991) 2181–2195. T.S. Ro, G.M. Homsy, J. Non-Newt. Fluid Mech. 57 (1995) 203–225. P.C. Huzyak, K.W. Koelling, J. Non-Newt. Fluid Mech. 71 (1997) 73–88.

1408 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]

F. Kamısßlı / International Journal of Engineering Science 43 (2005) 1388–1408

F. Kamisli, M.E. Ryan, Chem. Eng. J. 75 (1999) 167–176. F. Kamisli, M.E. Ryan, Chem. Eng. Sci. 56 (2001) 4913–4928. P.G. Saffmanand, G.I. Taylor, Proc. Roy. Soc. London A 245 (1958) 312–329. E. Pitts, J. Fluid Mech. 97 (1980) 53–64. J.M. McLean, P.G. Saffman, J. Fluid Mech. 102 (1981) 455–469. P. Tabeling, G. Zocchi, A. Libchaber, J. Fluid Mech. 177 (1987) 67–82. T. Bauman, T. Sullivan, S. Middleman, Chem. Eng. Commun. 14 (1982) 35–46. R.B. Bird, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, vol. I, John Wiley, New York, 1987, pp. 296–355. A.H. Nayfeh, Perturbation Methods, John Wiley, New York, 1973. M. VanDyke, Perturbation Methods in Fluid Mechanics, Parabolic Press, California, 1975. F. Kamisli, Chem. Eng. Proc. 42 (2003) 569–581. R.L. Fosdick, K.R. Rajagopal, Proc. Roy. Soc. London A 339 (1980) 351–377.