Back
 JAMP  Vol.9 No.8 , August 2021
Streamlines in the Two-Dimensional Spreading of a Thin Fluid Film: Blowing and Suction Velocity Proportional to the Height
Abstract: The two-dimensional spreading under gravity of a thin fluid film with suction (fluid leak-off) or blowing (fluid injection) at the base is considered. The thin fluid film approximation is imposed. The height of the thin film satisfies a nonlinear diffusion equation with a source/sink term. The Lie point symmetries of the nonlinear diffusion equation are derived and exist, which provided the fluid velocity at the base, vn satisfies a first order linear partial differential equation. The general form has algebraic time dependence while a special case has exponential time dependence. The solution in which vn is proportional to the height of the thin film is studied. The width of the base always increases with time even for suction while the height decreases with time for sufficiently weak blowing. The streamlines of the fluid flow inside the thin film are plotted by first solving a cubic equation. For sufficiently weak blowing there is a dividing streamline, emanating from the stagnation point on the centre line which separates the fluid flow into two regions, a lower region consisting of rising fluid and dominated by fluid injection at the base and an upper region consisting of descending fluid and dominated by spreading due to gravity. For sufficiently strong blowing the lower region expands to completely fill the whole thin film.

1. Introduction

Thin fluid films occur in nature and they also play a significant role in many technological processes such as in coating applications and in the spreading of paints. There are important applications of thin fluid films to the lung [1] and the eye [2]. Reviews of the research up to the year 2000 on the spreading of thin fluid films and liquid drops have been given by Oron et al. [3] and by Davis [4].

Thin fluid flows usually occur under the action of gravity [5], surface tension gradients [4] or on a rotating substrate [6] [7]. The moving contact line at the boundary of the thin fluid film has been studied by O’Brien and Schwartz [8] and by Hocking [9]. Davis and Hocking [10] [11] have considered the spreading of a liquid drop on a porous base. The spreading of an axisymmetric liquid drop with fluid injection or suction at the base has been investigated by Mason and Momoniat [12] who later extended their work to include surface tension [13].

A study of the streamlines of the flow in a thin film can illustrate how the thin fluid film spreads. Momoniat et al. [14] investigated the effect of slot injection or suction on the spreading of a thin film where the surface tension gradient and gravity effects are included. They also investigated the behaviour of the streamlines and how they evolve with suction or blowing. They found that a breakup of the streamlines occurs. Mason and Chung [15] considered the spreading of a thin liquid drop with slip at the base but without fluid injection or suction at the base. The streamlines made clear the properties of the flow. It was found that vortices can occur at the base and close to the axis of the liquid drop.

The motivation for this research is to obtain an improvement on our understanding of the fluid flow inside a spreading thin fluid film with suction or blowing at the base by analysing the streamlines of the flow.

In general a nonlinear diffusion equation for the height of the thin film occurs in the mathematical formulation of thin fluid film flows. Exact analytical solutions for these equations have been found in the form of similarity solutions [16] [17]. The application of Lie group analysis of differential equations has been successful in solving problems in thin fluid film theory [18] [19] [20] and this method will be used here to reduce the nonlinear diffusion equation obtained to an ordinary differential equation. The half-width of the base is obtained during this reduction.

In this paper we investigate the effect of suction and blowing at the base on the spreading under gravity of a two-dimensional thin fluid film. The general form of the suction/blowing velocity at the base is determined from the condition that the partial differential equation and boundary conditions admit an invariant solution in a Lie symmetry analysis. We consider the case in which the suction/blowing velocity at the base is proportional to the height of the thin fluid film. The streamlines of the fluid flow inside the thin film are derived and analysed.

A nomenclature is provided in Table 1.

Table 1. Nomenclature for significant variables and parameters.

2. Mathematical Model

Consider the two-dimensional spreading under gravity of a thin film of viscous incompressible fluid on a fixed horizontal plane. At the base there is either suction/leak-off of fluid from the thin film or blowing/injection of fluid into the thin film. The model is illustrated in Figure 1. The derivation of the thin fluid film equations is well documented in the literature [3] [8]. We follow the derivation of Acheson [21].

The thin film is infinite in the y-direction and symmetric in the ( y , z ) -plane. The equation of the free surface is

z = h ( t , x ) (2.1)

where

h ( t , x ) = h ( t , x ) (2.2)

and the width of the thin film is

w ( t ) x w ( t ) . (2.3)

The flow is two-dimensional in the ( x , z ) -plane with velocity components and pressure

v x = v x ( t , x , z ) , v y = 0 , v z = v z ( t , x , z ) , p = p ( t , x , z ) (2.4)

Figure 1. Two-dimensional thin fluid film showing coordinate system, body force and boundary conditions.

where

v x ( t , x , z ) = v x ( t , x , z ) , v z ( t , x , z ) = v z ( t , x , z ) , p ( t , x , z ) = p ( t , x , z ) . (2.5)

The x and z components of the Navier-Stokes equation and the conservation of mass equation are [21]

ρ ( v x t + v x v x x + v z v x z ) = p x + μ ( 2 v x x 2 + 2 v x z 2 ) , (2.6)

ρ ( v z t + v x v z x + v z v z z ) = p z ρ g + μ ( 2 v z x 2 + 2 v z z 2 ) , (2.7)

v x x + v z z = 0, (2.8)

where g is the acceleration due to gravity.

We introduce characteristic quantities and dimensionless variables and write (2.6) to (2.8) in dimensionless form:

characteristic length in the x-direction L = w ( 0 ) = w 0 ,

characteristic length in the z-direction H = h ( 0,0 ) ,

characteristic velocity in the x-direction U,

characteristic velocity in the z-direction V = H L U ,

characteristic pressure and stress P,

characteristic half-width w ( 0 ) = w 0 .

The characteristic velocity V is determined from the conservation of mass Equation (2.8) which is not approximated so that the two terms balance while U and P are yet to be specified. The Reynolds number is

Re = U L ν (2.9)

where ν = μ / ρ is the kinematic viscosity. Dimensionless variables are defined by

x ¯ = x L , z ¯ = z H , v ¯ x = v x U , v ¯ z = v z V , t ¯ = t T , p ¯ = p P , h ¯ = h H , w ¯ = w w 0 . (2.10)

The x-component of the Navier-Stokes Equation (2.6), becomes

Re ( H L ) 2 [ v ¯ x t ¯ + v ¯ x v ¯ x x ¯ + v ¯ z v ¯ x z ¯ ] = P H 2 μ L U p ¯ x ¯ + ( H L ) 2 2 v ¯ x x ¯ 2 + 2 v ¯ x z ¯ 2 . (2.11)

For the thin fluid film to spread in the x-direction the magnitude of the pressure gradient in the x-direction should be sufficiently large to balance the viscous term. Thus

P = μ L U H 2 . (2.12)

The z-component of the Navier-Stokes Equation (2.7), in dimensionless form is

Re ( H L ) 4 [ v ¯ z t ¯ + v ¯ x v ¯ z x ¯ + v ¯ z v ¯ z z ¯ ] = p ¯ z ¯ ρ g H P + ( H L ) 2 [ ( H L ) 2 2 v ¯ z x ¯ 2 + 2 v ¯ z z ¯ 2 ] , (2.13)

where (2.12) was used for P in all terms except the body force term. For the thin fluid film to spread the body force term must be sufficiently large to balance the pressure gradient term in the z-direction. Thus we obtain the second condition on P,

P = ρ g H . (2.14)

Equating (2.12) and (2.14) for P gives for the characteristic velocity U,

U = g H 3 ν L . (2.15)

The thin fluid approximation is [21]

Re ( H L ) 2 1 , ( H L ) 2 1. (2.16)

Equations (2.11) and (2.13) reduce to

2 v ¯ x z ¯ 2 = p ¯ x ¯ , (2.17)

p ¯ z ¯ = 1, (2.18)

and (2.8) expressed in dimensionless variables is

v ¯ x x ¯ + v ¯ z z ¯ = 0. (2.19)

Before stating the boundary conditions we derive two preliminary results. From the Navier-Poisson law for an incompressible fluid [21]

τ z z = p + 2 μ v z z , τ z x = μ ( v z x + v x z ) , (2.20)

where τ z z and τ z x are components of the Cauchy stress tensor. Expressed in dimensionless variables,

τ ¯ z z = p ¯ + 2 ( H L ) 2 v ¯ z z ¯ , τ ¯ z x = H L v ¯ x z ¯ + ( H L ) 2 v ¯ z x ¯ (2.21)

and neglecting terms of order ( H L ) 2 ,

τ ¯ z z = p ¯ , τ ¯ z x = H L v ¯ x z ¯ . (2.22)

Also, a fluid particle on the free surface of the thin film remains on the free surface as the fluid evolves. Thus

D D t ( z h ( t , x ) ) | z = h ( t , x ) = 0 (2.23)

where D D t is the convective time derivative. Hence

v z ( t , x , h ) = h t + v x ( t , x , h ) h x (2.24)

and expressed in terms of dimensionless variables,

v ¯ z ( t ¯ , x ¯ , h ¯ ) = h ¯ t ¯ + v ¯ x ( t ¯ , x ¯ , h ¯ ) h ¯ x ¯ . (2.25)

We now state the boundary conditions and then comment on them briefly.

z ¯ = 0 : v ¯ x ( t ¯ , x ¯ , 0 ) = 0 , (2.26)

z ¯ = 0 : v ¯ z ( t ¯ , x ¯ , 0 ) = v ¯ n ( t ¯ , x ¯ ) , (2.27)

z ¯ = h ¯ : p ¯ ( t ¯ , x ¯ , h ¯ ) = p ¯ 0 , (2.28)

z ¯ = h ¯ : τ ¯ z x ( t ¯ , x ¯ , h ¯ ) = H L v ¯ x z ¯ ( t ¯ , x ¯ , h ¯ ) = 0 , (2.29)

z ¯ = h ¯ : h ¯ t ¯ + v ¯ x ( t ¯ , x ¯ , h ¯ ) h ¯ x ¯ ( t ¯ , x ¯ ) = v ¯ z ( t ¯ , x ¯ , h ¯ ) . (2.30)

Boundary condition (2.26) is the no slip condition for a viscous fluid at a solid boundary. Boundary condition (2.27) balances the normal velocity induced in the thin fluid film with the velocity of fluid injection or suction at the base where

v ¯ n ( t ¯ , x ¯ ) > 0 for injection/blowing ,

v ¯ n ( t ¯ , x ¯ ) = 0 for non-porous base ,

v ¯ n ( t ¯ , x ¯ ) < 0 for suction .

The characteristic value for v n is the characteristic velocity in the z-direction, V, and

v n ( t , x ) = v n ( t , x ) . (2.31)

The boundary condition (2.28) is the balance of the normal stress at the surface of the thin film where p 0 is the atmospheric pressure and the boundary condition (2.29) is that the tangential stress vanishes at the surface.

Consider now the partial differential equation for the surface h ¯ ( t ¯ , x ¯ ) which is derived from (2.30). In order to obtain v ¯ z ( t ¯ , x ¯ , h ¯ ) , we integrate the conservation of mass Equation (2.19) with respect to z ¯ from z ¯ = 0 to z = h ¯ and use the boundary condition (2.27). This gives

v ¯ z ( t ¯ , x ¯ , h ¯ ) = v ¯ n ( t ¯ , x ¯ ) 0 h ¯ ( t ¯ , x ¯ ) v ¯ x x ¯ ( t ¯ , x ¯ , z ¯ ) d z ¯ . (2.32)

But using the formula for differentiation under the integral sign [22],

x ¯ 0 h ¯ ( t ¯ , x ¯ ) v ¯ x ( t ¯ , x ¯ , z ¯ ) d z ¯ = 0 h ¯ ( t ¯ , x ¯ ) v ¯ x x ¯ ( t ¯ , x ¯ , z ¯ ) d z ¯ + v ¯ x ( t ¯ , x ¯ , h ¯ ) h ¯ x ¯ (2.33)

and therefore (2.32) becomes

v ¯ z ( t ¯ , x ¯ , h ¯ ) = v ¯ n ( t ¯ , x ¯ ) + v ¯ x ( t ¯ , x ¯ , h ¯ ) h ¯ x ¯ x ¯ 0 h ¯ ( t ¯ , x ¯ ) v ¯ x ( t ¯ , x ¯ , z ¯ ) d z ¯ . (2.34)

Substituting (2.34) into (2.30) we obtain

h ¯ t ¯ + x ¯ 0 h ¯ ( t ¯ , x ¯ ) v ¯ x ( t ¯ , x ¯ , z ¯ ) d z ¯ = v ¯ n ( t ¯ , x ¯ ) . (2.35)

It remains to obtain v ¯ x ( t ¯ , x ¯ , z ¯ ) . Integrating (2.18) and imposing the boundary condition (2.28) gives

p ( t ¯ , x ¯ , z ¯ ) = z ¯ + h ¯ ( t ¯ , x ¯ ) + p ¯ 0 (2.36)

and Equation (2.17) becomes

2 v ¯ x z ¯ 2 = h ¯ x ¯ ( t ¯ , x ¯ ) . (2.37)

Integrating (2.37) twice with respect to z ¯ and imposing the no slip boundary condition (2.26) and the no tangential stress boundary condition (2.29) we obtain

v ¯ x ( t ¯ , x ¯ , z ¯ ) = ( z ¯ h ¯ ( t ¯ , x ¯ ) z ¯ 2 2 ) h ¯ x ¯ . (2.38)

Finally substituting (2.38) into (2.35) and integrating we find that

h ¯ t ¯ = 1 3 x ¯ ( h ¯ 3 h ¯ x ¯ ) + v ¯ n ( t ¯ , x ¯ ) . (2.39)

Equation (2.39) is a nonlinear diffusion equation for h ¯ ( t ¯ , x ¯ ) with a source/sink term v ¯ n ( t ¯ , x ¯ ) . The velocity v ¯ n ( t ¯ , x ¯ ) is not specified but is determined from the condition that the problem admits an invariant solution.

Since the height of the thin fluid film vanishes at the moving contact lines, x ¯ = ± w ¯ ( t ¯ ) , it follows that

h ¯ ( t ¯ , ± w ¯ ( t ¯ ) ) = 0. (2.40)

The volume flux of fluid in the x-direction per unit breadth in the y-direction is

Q ( t , x ) = 0 h ( t , x ) v x ( t , x , z ) d z . (2.41)

The characteristic volume flux is UH. Expressed in dimensionless variables and using (2.38),

Q ¯ ( t ¯ , x ¯ ) = 0 h ¯ ( t ¯ , x ¯ ) v ¯ x ( t ¯ , x ¯ , z ¯ ) = 1 3 h ¯ 3 ( t ¯ , x ¯ ) h ¯ x ¯ ( t ¯ , x ¯ ) . (2.42)

Consider now the balance law for fluid volume. The total volume of the thin fluid film per unit length in the y-direction is

V ¯ ( t ¯ ) = w ¯ ( t ¯ ) w ¯ ( t ¯ ) h ¯ ( t ¯ , x ¯ ) d x ¯ (2.43)

where the characteristic volume is H w 0 . Applying the formula for differentiation under the integral sign [22]

d V ¯ d t ¯ = w ¯ ( t ) w ¯ ( t ) h ¯ t ¯ ( t ¯ , x ¯ ) d x ¯ + h ¯ ( t ¯ , w ¯ ( t ) ) d w ¯ d t ¯ + h ¯ ( t ¯ , w ¯ ( t ¯ ) ) d w ¯ d t ¯ . (2.44)

Imposing the boundary conditions (2.40) and substituting the nonlinear diffusion Equation (2.39), (2.44) becomes

d V ¯ d t ¯ = Q ¯ ( t ¯ , w ¯ ( t ¯ ) ) Q ¯ ( t ¯ , w ¯ ( t ¯ ) ) + w ¯ ( t ) w ¯ ( t ) v ¯ n ( t ¯ , x ¯ ) d x ¯ . (2.45)

At the moving contact lines, x ¯ = ± w ¯ ( t ¯ ) , we take as boundary condition that the volume flux of fluid in the x-direction vanishes, that is

Q ¯ ( t ¯ , ± w ¯ ( t ¯ ) ) = 1 3 h ¯ 3 ( t ¯ , ± w ¯ ( t ¯ ) ) h ¯ x ¯ ( t ¯ , ± w ¯ ( t ¯ ) ) = 0. (2.46)

The balance law for fluid volume reduces to

d V d t = 2 0 w ¯ ( t ¯ ) v ¯ n ( t ¯ , x ¯ ) d x ¯ . (2.47)

The problem is to solve the nonlinear diffusion Equation (2.39) for h ¯ ( t ¯ , x ¯ ) subject to the boundary conditions (2.40) and (2.46) and to the requirement that (2.39) admits a Lie point symmetry which generates an invariant solution.

In the remainder of the paper the overhead bars are suppressed to keep the notation simple, it being understood that dimensionless quantities are used.

3. Lie Point Symmetries and Group Invariant Solutions

When deriving the Lie point symmetries of the nonlinear diffusion Equation (2.39) the variables t , x , h and all the partial derivatives of h with respect to t and x are regarded as independent variables. A suffix is used to denote the partial derivatives of h ( t , x ) . Equation (2.39) can be written in the form

F ( t , x , h , h t , h x , h x x ) = 0 (3.1)

where

F ( t , x , h , h t , h x , h x x ) = h t h 2 h x 2 1 3 h 3 h x x v n ( t , x ) . (3.2)

The Lie point symmetry generators of (2.39) are of the form [23]

X = ξ 1 ( t , x , h ) t + ξ 2 ( t , x , h ) x + η ( t , x , h ) h . (3.3)

They are derived from the determining equation

X [ 2 ] F | F = 0 = 0 (3.4)

where the second prolongation X [ 2 ] of X is defined by

X [ 2 ] = X + ζ 1 h t + ζ 2 h x + ζ 11 h t t + ζ 12 h t x + ζ 22 h x x (3.5)

where

ζ i = D i ( η ) h k D i ( ξ k ) , i = 1,2 (3.6)

ζ i j = D j ( ζ i ) h i k D j ( ξ k ) , i = 1,2 (3.7)

with summation over the repeated index k from 1 to 2. The total derivative operators with respect to t and x are

D 1 = D t = t + h t h + h t t h t + h x t h x + , (3.8)

D 2 = D x = x + h x h + h t x h t + h x x h x + . (3.9)

The determining Equation (3.4) is separated by powers and products of the partial derivatives of h. It is found that

X = ( c 1 + c 2 t ) t + ( c 3 + c 4 x ) x + 1 3 ( 2 c 4 c 2 ) h h = c 1 X 1 + c 2 X 2 + c 3 X 3 + c 4 X 4 (3.10)

where

X 1 = t , X 2 = t t 1 3 h h , X 3 = x , X 4 = x x 2 3 h h , (3.11)

provided the velocity v n ( t , x ) satisfies the first order linear partial differential equation

( c 1 + c 2 t ) v n t + ( c 3 + c 4 x ) v n x = 2 3 ( c 4 2 c 2 ) v n . (3.12)

The velocity v n ( t , x ) therefore cannot be arbitrary for the Lie point symmetries of Equation (2.39) to exist.

Now, h = Φ ( t , x ) is a group invariant solution of the nonlinear diffusion Equation (2.39) provided

X ( h Φ ( t , x ) ) | h = Φ ( t , x ) = 0 , (3.13)

that is, provided Φ ( t , x ) satisfies the first order linear partial differential equation

( c 1 + c 2 t ) Φ t + ( c 3 + c 4 x ) Φ x = 1 3 ( 2 c 4 c 2 ) Φ . (3.14)

We first consider the general case in which c 2 0 , c 4 0 and 2 c 4 c 2 0 and then consider the special case c 2 = 0 , c 4 0 and c 1 0 .

3.1. General Case c 2 0 , c 4 0 , c 2 2 c 4

The differential equations of the characteristic curves of (3.14) are

d t c 2 ( c 1 c 2 + t ) = d x c 4 ( c 3 c 4 + x ) = 3 d Φ ( 2 c 4 c 2 ) Φ . (3.15)

Integrating of the first pair of terms in (3.15) gives

c 3 c 4 + x ( c 1 c 2 + t ) c 4 / c 2 = a 1 , (3.16)

where a 1 is a constant while integration of the first and last pair of terms gives

Φ ( c 1 c 2 + t ) 2 3 ( c 4 c 2 1 2 ) = a 2 , (3.17)

where a 2 is a constant. The general solution of (3.14) is

a 2 = F ( a 1 ) (3.18)

where F is an arbitrary function. Hence, since h = Φ , the form of the invariant solution is

h ( t , x ) = ( c 1 c 2 + t ) 2 3 ( c 4 c 2 1 2 ) F ( ξ ) (3.19)

where

ξ = c 3 c 4 + x ( c 1 c 2 + t ) c 4 / c 2 . (3.20)

Consider next the partial differential Equation (3.12) for v n ( t , x ) . We assume in addition that c 4 2 c 2 . The differential equations of the characteristic curves of (3.12) are

d t c 2 ( c 1 c 2 + t ) = d x c 4 ( c 3 c 4 + x ) = 3 d v n 2 ( c 4 2 c 2 ) v n . (3.21)

Integration of the first pair of terms gives again (3.16) and integrating the first and last pair of terms gives

v n ( c 1 c 2 + t ) 2 3 ( c 4 c 2 2 ) = a 3 (3.22)

where a 3 is a constant. The general solution of (3.12) is

a 3 = G ( a 1 ) (3.23)

where G is an arbitrary function. Hence the form of the invariant solution is

v n ( t , x ) = ( c 1 c 2 + t ) 2 3 ( c 4 c 2 2 ) G ( ξ ) , (3.24)

where ξ is given by (3.20).

Substituting (3.19) and (3.24) into the nonlinear diffusion Equation (2.39) reduces (2.39) to the ordinary differential equation

1 3 d d ξ [ F 3 d F d ξ ] + c 4 c 2 d d ξ ( ξ F ) + 5 3 ( 1 5 c 4 c 2 ) F + G = 0. (3.25)

Equation (3.25) depends only on the ratio c 4 c 2 and does not depend on c 1 and c 3 .

The half-width of the base, w ( t ) , is determined from the boundary conditions in (2.40). Using (3.19) for h ( t , x ) , (2.40) becomes

F ( A ( t ) ) = 0 and F ( B ( t ) ) = 0 , (3.26)

where

A ( t ) = c 3 c 4 + w ( t ) ( c 1 c 2 + t ) c 4 / c 2 , B ( t ) = c 3 c 4 w ( t ) ( c 1 c 2 + t ) c 4 / c 2 . (3.27)

Differentiating (3.26) with respect to t gives

d F d A d A d t = 0 and d F d B d B d t = 0 (3.28)

and therefore since F ( ξ ) is not a constant function

d A d t = 0 and d B d t = 0. (3.29)

Thus A ( t ) = A 0 and B ( t ) = B 0 where A 0 and B 0 are constants and

c 3 c 4 + w ( t ) ( c 1 c 2 + t ) c 4 / c 2 = A 0 , c 3 c 4 w ( t ) ( c 1 c 2 + t ) c 4 / c 2 = B 0 . (3.30)

Adding we obtain

2 c 3 c 4 ( c 1 c 3 + t ) c 4 c 2 = A 0 + B 0 (3.31)

and differentiating (3.31) with respect to t, we find that

c 3 c 2 ( c 1 c 3 + t ) c 4 c 2 1 = 0. (3.32)

Thus c 3 = 0 and from (3.31), B 0 = A 0 . Equation (3.30) becomes

w ( t ) = A 0 ( c 1 c 2 + t ) c 4 c 2 (3.33)

and therefore since from (2.10), w ( 0 ) = 1 ,

A 0 = ( c 2 c 1 ) c 4 c 2 and w ( t ) = ( 1 + c 2 c 1 t ) c 4 c 2 . (3.34)

The boundary condition (3.26) becomes

F ( ± ( c 2 c 1 ) c 4 c 2 ) = 0 (3.35)

and from (3.20) with c 3 = 0 and (3.33),

ξ = ( c 2 c 1 ) c 4 c 2 x w ( t ) . (3.36)

The total volume of the thin fluid film per unit length in the y-direction is

V ( t ) = 2 0 w ( t ) h ( t , x ) d x . (3.37)

By using (3.19) and (3.36) it can be expressed as

V ( t ) = V 0 ( 1 + c 2 c 1 t ) 5 3 ( c 4 c 2 1 5 ) (3.38)

where

V 0 = 2 ( c 1 c 2 ) 5 3 ( c 4 c 2 1 5 ) 0 ( c 2 c 1 ) c 4 / c 2 F ( ξ ) d ξ . (3.39)

The balance law for fluid volume is given by (2.47) and using (3.24) for v n ( t , x ) it becomes

5 3 ( c 4 c 2 1 5 ) 0 ( c 2 c 1 ) c 4 / c 2 F ( ξ ) d ξ = 0 ( c 2 c 1 ) c 4 / c 2 G ( ξ ) d ξ . (3.40)

In order to simplify the ordinary differential Equation (3.25) and the boundary conditions (3.35) we make the change of variables

ξ = A η , F ( ξ ) = B f ( η ) , G ( ξ ) = D g ( η ) , (3.41)

where

A = ( c 2 c 1 ) c 4 c 2 , B = ( c 4 c 2 A 2 ) 1 3 , D = c 4 c 2 B . (3.42)

Equations (3.25) and (3.35) become

d d η ( f 3 d f d η ) + 3 d d η ( η f ) + ( c 2 c 4 5 ) f + 3 g = 0 (3.43)

and

f ( ± 1 ) = 0. (3.44)

The other variables transform to

h ( t , x ) = ( c 4 c 1 ) 1 3 ( 1 + c 2 c 1 t ) 2 3 ( c 4 c 2 1 2 ) f ( η ) , (3.45)

v n ( t , x ) = ( c 4 c 1 ) 4 3 ( 1 + c 2 c 1 t ) 2 3 ( c 4 c 2 2 ) g ( η ) , (3.46)

1 3 ( 5 c 2 c 4 ) 0 1 f ( η ) d η = 0 1 g ( η ) d η , (3.47)

V 0 = 2 ( c 4 c 1 ) 1 3 0 1 f ( η ) d η , (3.48)

where

η = x w ( t ) . (3.49)

There is one condition still to be imposed. Since the characteristic length in the z-direction is the initial height of the thin film at the central line

h ( 0,0 ) = 1. (3.50)

Hence from (3.45),

c 4 c 1 = 1 f 3 ( 0 ) . (3.51)

The ratio c 4 c 2 occurs in the differential Equation (3.43) and balance law (3.47) and is either given or determined from the balance law. Let

c 4 c 2 = α , (3.52)

then

c 2 c 1 = c 2 c 4 c 4 c 1 = 1 α f 3 ( 0 ) . (3.53)

The problem expressed in terms of the transformed variable η and the parameter α is to solve the differential Equation (3.43),

d d η ( f 3 d f η ) + 3 d d η ( η f ) + ( 1 α 5 ) f + 3 g = 0 (3.54)

subject to the boundary conditions (3.44) and (2.46),

f ( 1 ) = 0 , f ( 1 ) = 0 , (3.55)

f 3 ( 1 ) d f d η ( 1 ) = 0, f 3 ( 1 ) d f d η ( 1 ) = 0, (3.56)

and the balance law

1 3 ( 5 1 α ) 0 1 f ( η ) d η = 0 1 g ( η ) d η . (3.57)

The remaining quantities are given by

h ( t , x ) = ( 1 + t α f 3 ( 0 ) ) 2 3 ( α 1 2 ) f ( η ) f ( 0 ) , (3.58)

v n ( t , x ) = 1 f 4 ( 0 ) ( 1 + t α f 3 ( 0 ) ) 2 3 ( α 2 ) g ( η ) , (3.59)

w ( t ) = ( 1 + t α f 3 ( 0 ) ) α , (3.60)

V ( t ) = V 0 ( 1 + t α f 3 ( 0 ) ) 5 3 ( α 1 5 ) , (3.61)

where

V 0 = 2 0 1 f ( η ) f ( 0 ) d η (3.62)

and η is given by (3.49). The Lie point symmetry (3.10) which generates the invariant solution is

X = ( f 3 ( 0 ) + 1 α t ) t + x x + 1 3 ( 2 1 α ) h h . (3.63)

3.2. Special Case c 2 = 0 , c 4 0 and c 1 0

When c 2 = 0 the Lie point symmetry (3.10) reduces to

X = c 1 t + ( c 3 + c 4 x ) x + 2 3 c 4 h h . (3.64)

Now h = Φ ( t , x ) is a group invariant solution of (2.39) provided (3.13) is satisfied, that is, provided

c 1 Φ t + ( c 3 + c 4 x ) Φ x = 2 3 c 4 Φ . (3.65)

The differential equations of the characteristic curves of (3.65) are

d t c 1 = d x c 4 ( c 3 c 4 + x ) = 3 d Φ 2 c 4 Φ . (3.66)

The first pair of terms gives

c 3 c 4 + x exp ( c 4 c 1 t ) = b 1 (3.67)

while the pair consisting of the first and last terms give

Φ exp ( 2 3 c 4 c 1 t ) = b 2 , (3.68)

where b 1 and b 2 are constants. The general solution of (3.65) is

b 2 = F ( b 1 ) (3.69)

where F is an arbitrary function. Hence since h = Φ ,

h ( t , x ) = exp ( 2 3 c 4 c 1 t ) F ( ξ ) , (3.70)

where

ξ = c 3 c 4 + x exp ( c 4 c 1 t ) . (3.71)

When c 2 = 0 , equation (3.12) reduces to

c 1 v n t + ( c 3 + c 4 x ) v n x = 2 3 c 4 v n (3.72)

which has the same form as (3.65). Hence

v n ( t , x ) = exp ( 2 3 c 4 c 1 t ) G ( ξ ) (3.73)

where G ( ξ ) is an arbitrary function of ξ .

Substituting (3.70) and (3.73) into the partial differential Equation (2.39) reduces it to the ordinary differential equation

1 3 d d ξ ( F 3 d F d ξ ) + c 4 c 1 d d ξ ( ξ F ) 5 3 c 4 c 1 F + G = 0. (3.74)

The half-width of the base, w ( t ) , is again obtained from the boundary conditions (2.40) which can be written as (3.26) with

A ( t ) = c 3 c 4 + w ( t ) exp ( c 4 c 1 t ) , B ( t ) = c 3 c 4 w ( t ) exp ( c 4 c 1 t ) . (3.75)

It can be verified as for the general case that A ( t ) = A 0 and B ( t ) = B 0 where A 0 and B 0 are constants with B 0 = A 0 and that c 3 = 0 . Since w ( 0 ) = 1 we find that

w ( t ) = exp ( c 4 c 1 t ) . (3.76)

The similarity variable (3.71) becomes

ξ = x exp ( c 4 c 1 t ) = x w ( t ) (3.77)

and the boundary conditions (3.26) are

F ( 1 ) = 0 , F ( 1 ) = 0. (3.78)

The total volume (3.37), expressed in terms of the invariant solution, is

V ( t ) = V 0 exp ( 5 3 c 4 c 1 t ) , (3.79)

where

V 0 = 2 0 1 F ( ξ ) d ξ . (3.80)

The balance law (2.47) becomes

5 3 c 4 c 1 0 1 F ( ξ ) d ξ = 0 1 G ( ξ ) d ξ . (3.81)

The remaining condition, h ( 0,0 ) = 1 , yields

F ( 0 ) = 1. (3.82)

In order to treat the special case c 2 = 0 as the limit c 2 0 of the general case we transform the special case by making the change of variables (3.41) where now

A = 1 , B = ( c 4 c 1 ) 1 3 , D = ( c 4 c 1 ) 4 3 . (3.83)

Expressed in terms of the transformed variables the problem is to solve the differential equation

d d η ( f 3 d f d η ) + 3 d d η ( η f ) 5 f + 3 g = 0 (3.84)

subject to the boundary conditions (3.55) and (3.56) and to the balance law

5 3 0 1 f ( η ) d η = 0 1 g ( η ) d η (3.85)

and to the scaling condition

c 4 c 1 = 1 f 3 ( 0 ) . (3.86)

The remaining quantities are

h ( t , x ) = exp ( 2 3 t f 3 ( 0 ) ) f ( η ) f ( 0 ) , (3.87)

v n ( t , x ) = 1 f 4 ( 0 ) exp ( 2 3 t f 3 ( 0 ) ) g ( η ) , (3.88)

w ( t ) = exp ( t f 3 ( 0 ) ) , (3.89)

V ( t ) = V 0 exp ( 5 3 t f 3 ( 0 ) ) , (3.90)

where

V 0 = 2 0 1 f ( η ) f ( 0 ) d η , (3.91)

and

η = x w ( t ) . (3.92)

The Lie point symmetry which generates the invariant solution is

X = f 3 ( 0 ) t + x x + 2 3 h h . (3.93)

Unlike the Lie point symmetry (3.63) which generates the general solution, (3.93) is not a scaling symmetry.

The results (3.84) to (3.93) can be obtained from (3.54) to (3.63) by writing

( 1 + t α f 3 ( 0 ) ) α = exp [ α ln ( 1 + t α f 3 ( 0 ) ) ] (3.94)

and using the expansion

ln ( 1 + ε ) = ε + O ( ε 2 ) (3.95)

for | ε | < 1 we obtain

( 1 + t α f 3 ( 0 ) ) α = exp [ t f 3 ( 0 ) + O ( 1 α ) ] (3.96)

and therefore

lim 1 α 0 ( 1 + t α f 3 ( 0 ) ) α = exp ( t f 3 ( 0 ) ) . (3.97)

4. Normal Velocity at Base Proportional to the Thin Fluid Film Height

An assumption on v n ( t , x ) needs to be made to close the system of equations. In this section we will investigate the solution for which v n ( t , x ) is proportional to h ( t , x ) .

4.1. Invariant Solution

We will first consider the general case c 2 0 and then the special case c 2 = 0 . When c 2 0 suppose that

g ( η ) = β f ( η ) , < β < , (4.1)

where β is a constant. Then from (3.58) and (3.59)

v n ( t , x ) = β h ( t , x ) f 3 ( 0 ) ( 1 + t α f 3 ( 0 ) ) . (4.2)

Thus v n ( t , x ) is proportional to h ( t , x ) and it is greatest at the centre line and vanishes at the moving contact lines, x = ± w ( t ) .

The balance law (3.57) becomes

1 3 ( 5 1 α 3 β ) 0 1 f ( η ) d η = 0 (4.3)

and therefore

α = 1 5 3 β , β 5 3 . (4.4)

The value β = 5 3 and its relation to the special case c 2 = 0 will be considered later.

The differential Equation (3.54) reduces to

d d η ( f 3 d f d η ) + 3 d d η ( η f ) = 0. (4.5)

The solution of (4.5) which satisfies the boundary conditions (3.55) and (3.56) is

f ( η ) = ( 9 2 ) 1 3 ( 1 η 2 ) 1 3 . (4.6)

Hence using (3.58) to (3.62),

h ( t , x ) = [ 1 + 2 9 ( 5 3 β ) t ] β 1 5 3 β ( 1 η 2 ) 1 3 , (4.7)

v n ( t , x ) = 2 9 β [ 1 + 2 9 ( 5 3 β ) t ] 2 ( 2 β 3 ) ( 5 3 β ) ( 1 η 2 ) 1 3 , (4.8)

w ( t ) = [ 1 + 2 9 ( 5 3 β ) t ] 1 5 3 β , (4.9)

V ( t ) = V 0 [ 1 + 2 9 ( 5 3 β ) t ] β 5 3 β , (4.10)

V 0 = 2 0 1 ( 1 η 2 ) 1 3 d η = 2 5 π Γ ( 1 3 ) Γ ( 5 6 ) , (4.11)

where Γ ( n ) is the Gamma function [22] and

η = x w ( t ) = x [ 1 + 2 9 ( 5 3 β ) t ] 1 5 3 β . (4.12)

The Lie point symmetry (3.63) which generates the invariant solution becomes

X = ( 9 2 + ( 5 3 β ) t ) t + x x + ( β 1 ) h h (4.13)

Consider now the special case c 2 = 0 and that (4.1) is satisfied. Then from (3.87) and (3.88),

v n ( t , x ) = β f 3 ( 0 ) h ( t , x ) (4.14)

and therefore v n ( t , x ) is proportional to h ( t , x ) . Equation (4.14) agrees with (4.2) in the limit 1 α 0 . From balance law (3.85),

( 5 3 β ) 0 1 f ( η ) d η = 0 (4.15)

and therefore β = 5 3 . The differential Equation (3.84) reduces to (4.5) and the solution which satisfies the boundary conditions (3.55) and (3.56) is (4.6). Hence from (3.87) to (3.92) we obtain

h ( t , x ) = exp ( 4 27 t ) ( 1 η 2 ) 1 3 , (4.16)

v n ( t , x ) = 10 27 exp ( 4 27 t ) ( 1 η 2 ) 1 3 , (4.17)

w ( t ) = exp ( 2 9 t ) , (4.18)

V ( t ) = V 0 exp ( 10 27 t ) , (4.19)

where V 0 is given by (4.11) and

η = x w ( t ) = x exp ( 2 9 t ) . (4.20)

The Lie point symmetry which generates the invariant solution, (3.93), becomes

X = 9 2 t + x x + 2 3 h h . (4.21)

In the same way as (4.17) was derived it can be shown that

lim β 5 3 [ 1 + 2 9 ( 5 3 β ) t ] 1 5 3 β = exp ( 2 9 t ) (4.22)

and it can be verified that (4.16) to (4.20) can be obtained from the general results (4.7) to (4.12) as β 5 3 .

4.2. Time Evolution of the Invariant Solution

We now investigate the evolution of the thin fluid film with time. The behaviour of the fluid variables depends on the value of the constant β where < β < .

We first consider the second condition in the thin fluid film approximation (2.16) and determine the range of values of β for which it is satisfied. Return to dimensional variables and define

ε ( t ) = h ( t , 0 ) w ( t ) = h ( 0 , 0 ) w ( 0 ) h ¯ ( t ¯ , 0 ) w ¯ ( t ¯ ) = h ( 0 , 0 ) w ( 0 ) ε ¯ ( t ¯ ) (4.23)

where

ε ¯ ( t ¯ ) = h ¯ ( t ¯ ,0 ) w ¯ ( t ¯ ) , ε ¯ ( 0 ) = 1. (4.24)

We assume that the thin fluid film approximation is satisfied initially so that

h ( 0 , 0 ) w ( 0 ) 1. (4.25)

Hence the thin fluid film condition will be satisfied as the fluid evolves if either ε ¯ ( t ¯ ) remains constant ( ε ¯ ( 0 ) = 1 ) or ε ¯ ( t ¯ ) 0 as t ¯ or as t ¯ t 1 where t 1 , is some limiting time. If ε ¯ ( t ¯ ) as t ¯ or as t ¯ t 1 then the thin fluid film approximation will break down.

Return to dimensional variables. From (4.7) and (4.9) for β 5 3 ,

ε ( t ) = [ 1 2 3 ( β 5 3 ) t ] 1 3 ( β 2 β 5 3 ) (4.26)

and from (4.16) and (4.18) for β = 5 3 ,

ε ( t ) = exp ( 2 27 t ) . (4.27)

For < β < 5 3 , ε ( t ) 0 algebraically as t , for β = 5 3 , ε ( t ) 0 exponentially as t while for 5 3 < β < 2 , ε ( t ) 0 in the finite time

t 1 = 3 2 ( β 5 3 ) . (4.28)

For β = 2 , ε ( t ) remains constant at ε = 1 while for 2 < β < , ε ( t ) in the finite time t 1 . The thin fluid film approximation is therefore not satisfied for 2 < β < and the solution is not valid. For < β < 0 there is fluid leak-off or suction at the base and the fluid film remains thin as it evolves. For 0 < β < there is fluid injection or blowing at the base and as long as the blowing is not too strong ( 0 < β 2 ) the thin fluid film condition is satisfied. In the following we will consider only the range < β 2 .

Consider next the evolution with time of the half-width of the base, w ( t ) , given by (4.9) for β 5 3 and by (4.18) for β = 5 3 . Now

d w d t = { 2 9 [ 1 2 3 ( β 5 3 ) t ] 3 β 4 5 3 β > 0 , β 5 3 2 9 exp ( 2 9 t ) > 0 , β = 5 3 (4.29)

and therefore d w d t > 0 for < β < . The base therefore always expands even for < β < 0 which describes suction of fluid from the thin film. For < β < 5 3 , w ( t ) algebraically as t , w ( t ) exponentially as t for β = 5 3 and for 5 3 < β 2 , w ( t ) in the finite time t 1 defined by (4.28).

The total volume of the thin film per unit length in the y-direction, V ( t ) , can be analysed in the same way. The evolution of ε ( t ) , w ( t ) and V ( t ) is summarised in Table 2.

The remaining two quantities, h ( t , x ) and v n ( t , x ) have additional values of β to analyse. Consider first h ( t , x ) which is given by (4.7) for β 5 3 and by (4.16) for β = 5 3 . For all < β 2 ,

h x ( t , x ) = 2 3 w ( t ) β 3 x ( 1 x 2 w 2 ( t ) ) 2 3 (4.30)

and

h x ( t , x ) as x ± w ( t ) . (4.31)

Hence in the vicinity of the moving contact lines, x = ± w ( t ) , the thin fluid film approximation is no longer satisfied. The inclusion of blowing or suction does not eliminate the singularity at x = ± w ( t ) . Also, for < β 2 ,

Table 2. Time evolution of ε ( t ) , w ( t ) and V ( t ) for < β < . The finite time t 1 is defined by (4.28). The thin fluid film approximation is not satisfied for 2 < β .

h t ( t ,0 ) = 2 9 ( β 1 ) w ( t ) 2 ( 2 β 3 ) . (4.32)

For β < 1 , h t ( t , 0 ) < 0 and therefore the height of the centre line of the thin fluid film decreases with time even when there is blowing at the base, 0 < β < 1 . The rate of decrease in the height due to spreading is greater than the rate of increase due to blowing. When β = 1 the two rates of change balance and when β > 1 the rate of increase due to blowing is greater than the rate of decrease due to spreading. For < β < 1 , h ( t ,0 ) 0 algebraically as t . For β = 1 , h ( t ,0 ) remains constant, while for 1 < β < 5 3 , h ( t ,0 ) algebraically as t . For β = 5 3 , h ( t ,0 ) exponentially and for 5 3 < β 2 , h ( t ,0 ) in the finite time t 1 given by (4.28). The time evolution of h ( t ,0 ) is summarised in Table 3.

Finally v n ( t ,0 ) which describes suction and blowing at the base, is given by (4.8) for β 5 3 and by (4.17) for β = 5 3 . Now for β 5 3 ,

v n t ( t , 0 ) = 4 ( 2 9 ) 5 3 β ( β 3 2 ) [ 1 2 3 ( β 5 3 ) t ] 7 β 11 5 3 β . (4.33)

For < β < 0 , v n ( t , 0 ) < 0 which describes leak-off or suction at the base. For this range of β , v n t ( t , 0 ) > 0 and hence the magnitude of the suction decreases as t increases and v n ( t ,0 ) 0 algebraically as t . When β = 0 the base is impermeable. When β > 0 , v n ( t , 0 ) > 0 and there is fluid injection or blowing of fluid into the thin film at the base. For 0 < β < 3 2 , v n t ( t , 0 ) < 0

Table 3. Time evolution of h ( t ,0 ) for < β 2 . The finite time t 1 is defined by (4.28).

so that the blowing becomes weaker as time increases and v n ( t ,0 ) 0 algebraically as t . For β = 3 2 , v n ( t ,0 ) is constant and the blowing is constant in time. For 3 2 < β 2 , v n t ( t , 0 ) > 0 and the blowing becomes stronger as time increases. For 3 2 < β < 5 3 , v n ( t ,0 ) algebraically as t , for β = 5 3 , v n ( t ,0 ) exponentially as t and for 5 3 < β 2 , v n ( t ,0 ) in the finite time t 1 defined by (4.28). The results for the time evolution of v n ( t ,0 ) are summarised in Table 4.

In Figure 2 the evolution of h ( t ,0 ) with time for a range of values of β is presented. Similar figures are obtained for the other fluid variables.

4.3. Stagnation Points on the Centre Line

In this section we investigate the z-component of the fluid velocity on the centre line, v z ( t ,0, z ) . Since from symmetry v x ( t ,0, z ) = 0 , points on the centre line at which v z ( t ,0, z ) = 0 are stagnation points.

To obtain v z ( t , x , z ) we integrate the continuity Equation (2.19) with respect to z from 0 to z and impose the boundary condition (2.27) which gives

Table 4. Time evolution of v n ( t ,0 ) for < β 2 . The finite time t 1 is defined by (4.28).

Figure 2. The maximum height of the thin fluid, h ( t ,0 ) , plotted against t for a range of values of β .

v z ( t , x , z ) = v n ( t , x ) x 0 z v x ( t , x , z ) d z . (4.34)

By using (2.38) for v x ( t , x , z ) we obtain

v z ( t , x , z ) = 1 6 z 3 2 h x 2 ( t , x ) + 1 4 z 2 2 x 2 ( h 2 ( t , x ) ) + v n ( t , x ) (4.35)

and therefore on the centre line x = 0 ,

v z ( t , 0 , z ) = 1 6 z 3 2 h x 2 ( t , x ) | x = 0 + 1 4 z 2 2 x 2 ( h 2 ( t , x ) ) | x = 0 + v n ( t , 0 ) . (4.36)

But from (4.7) and (4.9) for β 5 3 and from (4.16) and (4.18) for β = 5 3

h ( t , x ) = h ( t ,0 ) ( 1 x 2 w 2 ( t ) ) 1 3 (4.37)

where

h ( t , 0 ) = w ( t ) β 1 , h ( 0 , 0 ) = 1 (4.38)

and on the centre line, 0 z h ( t ,0 ) . Thus

2 h x 2 ( t , x ) | x = 0 = 2 3 h ( t , 0 ) w 2 ( t ) , (4.39)

2 x 2 ( h 2 ( t , x ) ) | x = 0 = 4 3 h 2 ( t , 0 ) w 2 ( t ) . (4.40)

It remains to obtain v n ( t ,0 ) . From (4.8) and (4.9) for β 5 3 and (4.17) and (4.18) for β = 5 3 ,

v n ( t , x ) = 2 9 β w ( t ) 2 ( 2 β 3 ) ( 1 η 2 ) 1 3 (4.41)

and therefore using (4.38)

v n ( t ,0 ) = 2 9 β h 4 ( t ,0 ) w 2 ( t ) . (4.42)

Substituting (4.39), (4.40) and (4.42) into (4.36) gives

v z ( t ,0, z ) = 1 9 h 4 ( t ,0 ) w 2 ( t ) [ z 3 3 z 2 + 2 β ] (4.43)

where

z * = z h ( t , 0 ) , 0 z * 1. (4.44)

The points on the centre line x = 0 at which v z ( t ,0, z ) = 0 satisfy the cubic equation

z 3 3 z 2 + 2 β = 0 , < β 2. (4.45)

The roots of the cubic Equation (4.45) which satisfy 0 z * 1 are stagnation points of the flow.

We first transform (4.45) to standard form [24]

s 3 + 3 H s + G = 0 (4.46)

by making the transformation

z * = s + 1. (4.47)

The cubic Equation (4.45) becomes

s 3 3 s 2 ( 1 β ) = 0. (4.48)

Comparison of (4.48) with (4.46) gives

H = 1 , G = 2 ( 1 β ) . (4.49)

The roots of (4.48) which satisfy 1 s 0 are stagnation points of the flow. The discriminant Δ of the cubic Equation (4.48) is

Δ = G 2 + 4 H 3 = 4 β ( β 2 ) . (4.50)

The roots of (4.48) satisfy the following general properties [24].

< β < 0 : Δ > 0 1 real root , 2 complex conjugate roots , β = 0 : Δ = 0 3 real roots , 2 coincident roots , 0 < β < 2 : Δ < 0 3 real and distinct roots, β = 2 : Δ = 0 3 real roots, 2 coincident roots, 2 < β < : Δ > 0 1 real root, 2 complex conjugate roots .

The thin fluid film approximation breaks down for 2 < β < .

The real zero for Δ > 0 does not lie in the range 1 s 0 . To show this let

P ( s ) = s 3 3 s 2 ( 1 β ) . (4.51)

Consider first < β < 0 . Then P ( 2 ) = 2 β < 0 and P ( ) = + . The real zero therefore lies in the range 2 < s < . Consider next 2 < β < . Then P ( 1 ) = 2 β > 0 and P ( ) = . The real zero therefore lies in range < s < 1 . We therefore consider the real roots of (4.48) in the range 0 β 2 for which Δ 0 . We use the trigonometric solution of the cubic equation when Δ 0 because the Cardan method gives the solution in complex form even although the roots are real [24].

We look for a solution of the form

s = ρ cos θ , ρ > 0. (4.52)

Substituting (4.52) into (4.48) gives

cos 3 θ + 3 H ρ 2 cos θ + G ρ 3 = 0. (4.53)

Comparing (4.53) with the trigonometric identity

cos 3 θ 3 4 cos θ 1 4 cos 3 θ = 0 (4.54)

we find that

ρ = 2 ( H ) 1 2 , cos 3 θ = G 2 ( H ) 3 2 . (4.55)

The condition 1 cos 3 θ 1 is satisfied because

G 2 + 4 H 3 = ρ 6 16 ( cos 3 3 θ 1 ) 0. (4.56)

Let ϕ be the smallest non-negative angle which satisfies

cos 3 ϕ = G 2 ( H ) 3 2 = 1 β (4.57)

and 0 3 ϕ π , that is, 0 ϕ π 3 . The general solution for θ is

θ = ϕ + 2 n π 3 , n = 0 , 1 , 2. (4.58)

Hence the three real roots of (4.48) for 0 β 2 are, from (4.52) and (4.55),

s n = 2 ( H ) 1 2 cos ( ϕ + 2 n π 3 ) , n = 0 , 1 , 2 (4.59)

and since z * = s + 1 and H = 1 the three real roots of the cubic Equation (4.45) are

z n * = 1 + 2 cos ( ϕ + 2 n π 3 ) , n = 0 , 1 , 2. (4.60)

4.3.1. Impermeable Base β = 0

From (4.57), cos 3 ϕ = 1 and ϕ = 0 since 0 ϕ π 3 . Thus

z 0 * = 3 , z 1 * = 0 , z 2 * = 0 , (4.61)

which may also be derived directly from (4.45) with β = 0 . The root z 0 * lies outside the thin fluid film. That (0,0) is a stagnation point follows directly from the boundary conditions.

4.3.2. Weak Blowing 0 < β < 1

We will derive expansions for z 0 * , z 1 * and z 2 * for small values of β . Since ϕ = 0 when β = 0 we expect ϕ to be small and positive when β is small. From (4.57)

cos 3 ϕ = 1 9 2 ϕ 2 + 27 8 ϕ 4 + O ( ϕ 6 ) = 1 β (4.62)

and neglecting terms O ( ϕ 6 ) we obtain form (4.62) the quadratic equation for ϕ 2 .

ϕ 4 4 3 ϕ 2 + 8 27 β = 0. (4.63)

Hence

ϕ = 2 3 β 1 2 ( 1 + O ( β ) ) . (4.64)

From (4.60) with n = 0 ,

z 0 * = 1 + 2 cos ϕ (4.65)

and

cos ϕ = 1 ϕ 2 2 + O ( ϕ 4 ) = 1 1 9 β + O ( β 2 ) . (4.66)

Thus

z 0 * = 3 ( 1 2 27 β + O ( β 2 ) ) (4.67)

which, although decreasing in magnitude lies outside the thin fluid film.

From (4.60) with n = 1 ,

z 1 * = 1 + 2 cos ( ϕ + 2 π 3 ) = 1 cos ϕ 3 sin ϕ . (4.68)

But

sin ϕ = ϕ 1 6 ϕ 3 + O ( ϕ 5 ) = 2 3 β 1 2 + O ( β 3 2 ) (4.69)

and substituting (4.66) and (4.69) into (4.68) we obtain

z 1 * = ( 2 3 ) 1 2 β 1 2 + 1 9 β + O ( β 3 2 ) . (4.70)

The root z 1 * is negative and lies outside the thin fluid film.

From (4.60) with n = 2 ,

z 2 * = 1 + 2 cos ( ϕ + 4 π 3 ) = 1 cos ϕ + 3 sin ϕ (4.71)

and substituting (4.66) and (4.69) into (4.71) we obtain

z 2 * = ( 2 3 ) 1 2 β 1 2 + 1 9 β + O ( β 3 2 ) . (4.72)

The root z 2 * lies inside the thin fluid film and the point ( 0, z 2 * ) is a stagnation point on the centre line x = 0 .

In Table 5 the approximate solution (4.72) is compared with the exact solution (4.71) where ϕ ( 0 ϕ π 3 ) satisfies (4.57) for 0 β 1 .

The expansions (4.67), (4.70) and (4.72) can also be derived using perturbation methods [25]. The solution (4.67) is a regular perturbation expansion while (4.70) and (4.72) are singular perturbation expansions.

4.3.3. Moderate Blowing β = 1

From (4.57) with β = 1 , cos 3 ϕ = 0 and since 0 ϕ π 3 it follows that ϕ = π 6 . Hence from (4.60)

z 0 * = 1 + 3 , z 1 * = 1 3 , z 2 * = 1. (4.73)

The roots z 0 * and z 1 * lie outside the thin fluid film. The root z 2 * = 1 is the limiting value of the root lying inside the thin film. For β > 1 , the root z 2 * > 1 and lies outside the thin film. The point ( 0, h ( t ,0 ) ) is the limiting case of a stagnation point on the centre line. The maximum height of the thin film remains constant with time. This is consistent with the results of Table 3.

4.3.4. Strong Blowing β = 2

When β = 2 , from (4.57), cos 3 ϕ = 1 and since 0 ϕ π 3 we have ϕ = π 3 . Thus from (4.60),

z 0 * = 2 , z 1 * = 1 , z 2 * = 2 , (4.74)

which can be derived directly from (4.45) with β = 2 . There are no roots inside the thin fluid film for 1 < β 2 .

The stagnation point z 2 * on the centre line plays a significant part in the pattern of the streamlines in the thin film which we now consider.

Table 5. The exact and approximate solution for the zero, z 2 / h ( t ,0 ) , of v z ( t ,0, z ) on the centre line.

4.4. Streamlines

In order to study further the fluid flow in the thin film we investigate the streamlines of the flow. The tangent vector to a streamline is everywhere parallel to the fluid velocity vector instantaneously. In two-dimensional incompressible flow the streamlines can be obtained from the stream function ψ ( t , x , z ) which satisfies

v x ( t , x , z ) = ψ z , v z ( t , x , z ) = ψ x . (4.75)

Since the stream function is constant along a streamline the streamlines in the ( x , z ) plane at time t can be obtained by plotting

ψ ( t , x , z ) = k ( t ) (4.76)

for a range of values of k ( t ) .

By substituting (2.38) and (4.36) for v x ( t , x , z ) and v z ( t , x , z ) into (4.75) we obtain

ψ z = ( 1 2 z 2 z h ( t , x ) ) h x ( t , x ) . (4.77)

ψ x = 1 6 z 3 2 h x 2 ( t , x ) 1 4 z 2 2 x 2 ( h 2 ( t , x ) ) v n ( t , x ) . (4.78)

It is readily verified that the compatibility condition

2 ψ x z = 2 ψ z x (4.79)

is satisfied. Integrating (4.77) yields

ψ ( t , x , z ) = 1 6 h x z 3 1 2 h h x z 2 + M ( t , x ) (4.80)

and by substituting (4.80) into (4.78) we find that

M x ( t , x ) = v n ( t , x ) . (4.81)

Hence

M ( t , x ) = 0 x v n ( t , s ) d s + D ( t ) (4.82)

where D ( t ) = M ( t ,0 ) . Thus

ψ ( t , x , z ) = A ( t , x ) z 3 + B ( t , x ) z 2 + C ( t , x ) + D ( t ) (4.83)

where

A ( t , x ) = 1 6 h x ( t , x ) , B ( t , x ) = 1 2 h ( t , x ) h x , C ( t , x ) = 0 x v n ( t , s ) d s . (4.84)

The streamlines in the ( x , z ) plane at time t are given by (4.76). They are obtained by solving numerically for z as a function of x the cubic equation

A ( t , x ) z 3 + B ( t , x ) z 2 + C ( t , x ) = K ( t ) (4.85)

at a given time t and a range of values of K ( t ) = k ( t ) D ( t ) . The cubic Equation (4.85) for the streamlines is quite general and does not depend on any relation between v n ( t , x ) and h ( t , x ) .

Consider now the invariant solutions (4.7) to (4.9) and (4.16) to (4.18). For both β 5 3 and β = 5 3 ,

A ( t , x ) = 1 9 w ( t ) β 3 x ( 1 x 2 w 2 ( t ) ) 2 3 , (4.86)

B ( t , x ) = 1 3 w ( t ) 2 ( β 2 ) x ( 1 x 2 w 2 ( t ) ) 1 3 , (4.87)

C ( t , x ) = 2 9 β w ( t ) 2 ( 2 β 3 ) 0 x ( 1 s 2 w 2 ( t ) ) 1 3 d s . (4.88)

The range of the base is w ( t ) x w ( t ) .

The streamlines derived from (4.85) depend on time. To determine the evolution of the streamline pattern with time it would have to be plotted for a range of values of time. Instead we will plot the streamlines in the ( η , z * ) plane where

η = x w ( t ) ( 1 η 1 ) , z * = z h ( t , 0 ) ( 0 z * 1 ) (4.89)

which we will see are independent of time and depend only on the parameter β . Clearly the base and maximum height are independent of time. The results for the evolution in time of w ( t ) and h ( t ,0 ) derived in Section 4.2 can then be used to determine how the streamline pattern evolves with time.

The cubic Equation (4.85) expressed in terms of η and z * is

A ( t , η ) z 3 + B ( t , η ) z 2 + C ( t , η ) = K ( t ) (4.90)

where using (4.86) to (4.88) and also (4.38)

A * ( t , η ) = h 3 ( t , 0 ) A ( t , x ) = 1 9 w ( t ) 4 β 5 η ( 1 η 2 ) 2 3 , (4.91)

B * ( t , η ) = h 2 ( t , 0 ) B ( t , x ) = 1 3 w ( t ) 4 β 5 η ( 1 η 2 ) 1 3 , (4.92)

C * ( t , η ) = C ( t , x ) = 2 9 β w ( t ) 4 β 5 0 η ( 1 ζ 2 ) 1 3 d ζ . (4.93)

Equations (4.91) to (4.93) are satisfied for both β 5 3 and β = 5 3 . The coefficients A * , B * and C * have the same time dependent factor w ( t ) 4 β 5 . Hence the cubic Equation (4.90) may be written as

η ( 1 η 2 ) 2 3 z 3 3 η ( 1 η 2 ) 1 3 z 2 + 2 β 0 η ( 1 ζ 2 ) 1 3 d ζ = K (4.94)

where

K * = 9 K ( t ) w ( t ) 4 β 5 . (4.95)

The streamlines are obtained by giving K * a range of values. Since (4.94) does not depend on time explictly the streamline pattern in the ( η , z * ) plane depends only on the parameter β and is independent of time.

In Figure 3 the streamlines are plotted in the ( η , z * ) plane for an impermeable base with β = 0 and for fluid leak-off (suction) at the base with β = 0.5 . For β = 0 the spreading is due only to gravity. The streamlines begin and end on the free surface except the streamline η = 0 which is the centre line. The point (0,0) is a stagnation point because of the boundary conditions. For fluid leak-off or suction the streamlines begin on the free surface and either end on the free surface or on the base. For stronger suction described by more negative values of β the streamlines are less curved and more streamlines end on the base. The streamlines are always perpendicular to the base because there is no slip at the base. There are no stagnation points in the flow.

In Figure 4 the streamlines are plotted in the ( η , z * ) plane for values of β in the range 0 < β 1 which describe fluid injection into the thin film at the base. There is a stagnation point on the centre line at ( 0, z 2 * ) where

z 2 * = 1 + 2 cos ( ϕ + 4 π 3 ) = ( 2 3 ) 1 2 β 1 2 + 1 9 β + O ( β 3 2 ) , (4.96)

where ϕ is given by (4.57). Dividing streamlines emanate from the stagnation point ( 0, z 2 * ) in the positive and negative η -directions and divide the flow into two regions. In the upper region the streamlines begin and end on the free surface except the streamline along the centre line which ends at the stagnation point. The flow is similar to that for β = 0 and is dominated by spreading due to gravity. In the lower region the flow is dominated by fluid injection at the base. This is especially the case near the centre line and for small values of z * but as z * increases spreading due to gravity causes the streamlines to bend towards the free surface. The maximum thickness of the lower region is z 2 * . As the strength of the blowing, β , increases we see that the thickness of the lower region increases. It forms a larger fraction of the thin film in agreement with the

Figure 3. Streamlines for (a) impermeable base with β = 0 and (b) leak-off at the base with β = 0.5 , plotted in the ( η , z * ) plane.

Figure 4. Streamlines for (a) β = 0.1 , (b) β = 0.2 , (c) β = 0.4 , (d) β = 0.6 , (e) β = 0.8 and (f) β = 1 , plotted in the ( η , z * ) plane.

numerical vales in Table 5. When β increases to β = 1 the lower region completely fills the thin film and there is no longer two regions. When β = 1 , the fluid velocity on the centre line is positive except at z * = 1 where v z ( t ,0,1 ) = 0 .

In Figure 5 the streamlines are plotted in the ( η , z * ) plane for β = 1.5 , β = 5 3 and for the limit of the thin fluid film approximation, β = 2 . The streamline pattern is similar to that for β = 1 . The streamlines emanate from the base at right angles because of the no slip boundary condition. Spreading due to gravity causes them to bend towards the free surface. The blowing is stronger

Figure 5. Streamlines for (a) β = 1.5 , (b) β = 5 3 and (c) β = 2 , plotted in the ( η , z * ) plane.

as β increases and the bending occurs at larger values of z * . On the centre line v z ( t ,0, z * ) > 0 for 0 z * 1 and there is no stagnation point in the flow.

The evolution of the streamline pattern with time can be determined using

x = w ( t ) η , z = h ( t ,0 ) z * (4.97)

and the results in Table 2 and Table 3 for the evolution of ε ( t ) , w ( t ) and h ( t ,0 ) with time.

In Figure 3, < β 0 and w ( t ) and h ( t ,0 ) 0 as t . The base will steadily increase and the maximum height will steadily decrease for 0 < t < .

In Figure 4, 0 < β 1 and w ( t ) as t . When 0 < β < 1 , h ( t ,0 ) 0 as t but for β = 1 , h ( t ,0 ) remains constant and the point of maximum height is a stagnation point. The height of the stagnation points above the base on the centre line,

z 2 = h ( t ,0 ) z 2 * (4.98)

will decrease with time for 0 β < 1 and remain constant for β = 1 . The maximum thickness of the lower and upper regions will both decrease with time but the ratio of the maximum thicknesses will remain constant as time increases:

Figure 6. Streamlines for β = 0.4 for (a) t = 1 and (b) t = 10 , plotted in the ( x , z ) -plane.

maximum thickness of upper region maximum thickness of lower region = h ( t ,0 ) z 2 z 2 = 1 z 2 * z 2 * . (4.99)

For β = 1 , the maximum height h ( t ,0 ) will remain constant but the base w ( t ) will expand to + as t .

In Figure 5, 1 < β 2 and w ( t ) and h ( t ,0 ) algebraically as t for 1 < β < 5 3 , exponentially for β = 5 3 and in the finite time t 1 defined by (4.28) for 5 3 < β 2 . The thin fluid film ratio ε ( t ) 0 as t or t t 1 for 1 < β < 2 and therefore the base expands to infinity faster than the maximum height but for β = 2 , ε ( t ) 1 and w ( t ) and h ( t ,0 ) tend to infinity at the same rate. For 2 < β < the thin fluid film approximation is not satisfied for all time.

The evolution of the streamlines with time is illustrated in Figure 6 where the streamlines are plotted in the ( x , z ) -plane for β = 0.4 at times t = 1 and t = 10 by solving numerically the cubic Equation (4.85) where A ( t , x ) , B ( t , x ) and C ( t , x ) are given by (4.86) to (4.88).

5. Conclusions

We first list the findings of the paper and then comment on the findings.

The thin fluid film equations were derived and it was shown how the characteristic fluid pressure could be determined from the dimensionless partial differential equations and how the characteristic horizontal fluid velocity could be determined from the two expressions obtained for the characteristic fluid pressure.

It was found that the height of the thin fluid film satisfied a nonlinear diffusion equation with a source/sink term consisting of the fluid injection/leak-off velocity, v n ( t , x ) , at the base.

The nonlinear diffusion equation for the height of the thin fluid film was found to possess a Lie point symmetry provided v n ( t , x ) which satisfies a first order linear partial differential equation. This leads to the invariant form for v n ( t , x ) .

The general form of the Lie point symmetry generates an invariant solution with height and half-width which evolve algebraically with time. A special form of the Lie point symmetry was also found which generates a solution for which the height and half-width evolve exponentially with time.

To close the system of equations, we found that one further condition is required. We assumed that v n ( t , x ) is proportional to the height h ( t , x ) .

The time evolution of the invariant solution for the thin fluid film ratio ε ( t ) , the height h ( t , x ) , the total volume per unit breadth V ( t ) and v n ( t , x ) was determined for both the general and special solution.

The stagnation point on the centre line was found by solving a cubic equation. At this point the upward flow due to fluid injection at the base is balanced by the downward flow due to spreading by gravity.

The fluid flow inside the thin fluid film was investigated by plotting the streamlines which were obtained from the stream function by solving numerically a second cubic equation. Streamlines for both fluid leak-off and fluid injection were plotted. For fluid injection, it was found that, provided the injection velocity is not too strong, there is a dividing streamline which separates the flow into two regions, an upper region with downflow due to gravity and a lower region with upflow due to fluid injection at the base.

We make the following comments on the findings.

When the fluid velocity at the base is proportional to the height of the thin fluid film, a complete analytical solution can be derived and the fluid variables and streamlines can be fully analysed. The thin fluid film approximation is satisfied for all values of time for the whole range of suction and significantly it is also satisfied for blowing at the base which provided the blowing that is not too strong.

There was unexpected behaviour of fluid variables due to the relative importance of spreading by gravity compared with blowing or suction at the base. The base half-width increased with time for all values of β , even for suction when < β < 0 . The maximum height of the thin film decreased for all time for suction and even for blowing provided the blowing was not too strong ( 0 < β < 1 ). Some limits were attained in a finite time and in such a way that the thin fluid film approximation remained satisfied. For sufficiently strong blowing ( 5 3 < β 2 ) the base half-width and the maximum height of the thin fluid film tended to infinity in time t 1 .

Much of the literature on thin fluid films is concerned with how the surface profile evolves with time. By plotting the streamlines numerically we were able to investigate the fluid flow inside the thin film. The streamline pattern for blowing was much richer than that for suction which did not show any unexpected features. The depth of penetration of the injected fluid at the base and its effect on the fluid flow inside the thin film could be investigated. The streamline pattern showed the existence of a dividing streamline which separated the flow into two regions, a lower region at the base consisting of rising fluid and an upper region consisting of descending fluid. We can expect that this structure will exist for other models of fluid injection. A streamline pattern independent of time was obtained. This was useful because from it and the known time evolution of w ( t ) and h ( t ,0 ) the streamline pattern at any later time could be determined.

The general solution for β 5 3 and the special solution for β = 5 3 where both generated by the Lie point symmetry (3.10). The general solution was generated from (3.10) with c 2 0 and the special solution with c 2 = 0 . This unified the two solutions and is the reason why the special solution could be obtained as a limiting case of the general solution. When c 2 0 , the Lie point symmetry is a scaling symmetry while the special case c 2 = 0 is not a scaling symmetry.

The theory developed in the paper up to the end of Section 3 is quite general. It is only in Section 4 that the assumption that v n ( t , x ) is proportional to h ( t , x ) is made. Other models for v n can be considered that lead to analytical solutions for the streamlines, for example, v n proportional to the spatial gradient of the height which was considered in the axisymmetric spreading of a liquid drop [12].

Acknowledgements

DPM thanks the National Research Foundation, Pretoria, South Africa, for financial support under grant number 132189.

Cite this paper: Modhien, N. , Mason, D. and Momoniat, E. (2021) Streamlines in the Two-Dimensional Spreading of a Thin Fluid Film: Blowing and Suction Velocity Proportional to the Height. Journal of Applied Mathematics and Physics, 9, 2114-2151. doi: 10.4236/jamp.2021.98133.
References

[1]   Grotberg, J.B. (1994) Pulmonary Flow and Transport Phenomena. Annual Review of Fluid Mechanics, 26, 529-571.
https://doi.org/10.1146/annurev.fl.26.010194.002525

[2]   Wong, H., Fatt, I. and Radke, C.J. (1996) Deposition and Thinning of the Human Tear Film. Journal of Colloid and Interface Science, 184, 44-51.
https://doi.org/10.1006/jcis.1996.0595

[3]   Oron, A., Davis, S.H. and Bankoff, S.G. (1997) Long-Scale Evolution of Thin Liquid Films. Reviews of Modern Physics, 69, 931-980.
https://doi.org/10.1103/RevModPhys.69.931

[4]   Davis, S.H. (2000) Interfacial Fluid Dynamics. In: Batchelor, G.K., Moffatt, M.K. and Worster, M.G., Eds., Perspectives in Fluid Dynamics, Cambridge University Press, Cambridge, 1-51.

[5]   Huppert, H.E. (1982) The Propagation of Two-Dimensional and Axisymmetric Viscous Gravity Currents over a Rigid Horizontal Surface. Journal of Fluid Mechanics, 121, 43-58.
https://doi.org/10.1017/S0022112082001797

[6]   Momoniat, E. and Mason, D.P. (1998) Investigation of the Effect of the Coriolis Force on a Thin Liquid Film on a Rotating Disk. International Journal of Non-Linear Mechanics, 33, 1069-1088.
https://doi.org/10.1016/S0020-7462(97)00071-1

[7]   Momoniat, E. and Myers, T.G. (2006) A New Solution for Rotation-Driven Spreading of a Thin Fluid Film. International Journal of Non-Linear Mechanics, 41, 102-109.
https://doi.org/10.1016/j.ijnonlinmec.2005.07.001

[8]   O’Brien, S.B.G. and Schwartz, L.W. (2006) Thin Film Flows: Theory and Modelling. In: Encyclopedia of Surface and Colloid Science, Taylor and Francis, Abingdon-on-Thames, 6304-6317.

[9]   Hocking, L.M. (1992) Rivel Contact-Angle Models and the Spreading of Drops. Journal of Fluid Mechanics, 239, 671-681.
https://doi.org/10.1017/S0022112092004579

[10]   Davis, S.H. and Hocking, L.M. (1999) Spreading and Imbibitions of Viscous Liquid on a Porous Base. Physics of Fluids, 11, 48-57.
https://doi.org/10.1063/1.869901

[11]   Davis, S.H. and Hocking, L.M. (2000) Spreading and Imbibitions of Viscous Liquid on a Porous Base II. Physics of Fluids, 12, 1646-1655.
https://doi.org/10.1063/1.870416

[12]   Mason, D.P. and Momoniat, E. (2004) Axisymmetric Spreading of a Thin Liquid Drop with Suction or Blowing at the Horizontal Base. International Journal of Non-Linear Mechanics, 39, 1013-1026.
https://doi.org/10.1016/S0020-7462(03)00093-3

[13]   Momoniat, E. and Mason, D.P. (2007) Spreading of a Thin Film with Suction or Blowing Including Surface Tension Effects. Computers & Mathematics with Applications, 53, 198-208.
https://doi.org/10.1016/j.camwa.2006.02.019

[14]   Momoniat, E., Ravindarin, R. and Roy, S. (2010) The Influence of Slot Injection/Suction on the Spreading of a Thin Film under Gravity and Surface Tension. Acta Mechanica, 211, 61-71.
https://doi.org/10.1007/s00707-009-0215-y

[15]   Mason, D.P. and Chung, M.T. (2007) Group Invariant Solution for Axisymmetric Spreading of a Thin Fluid Film with Slip at the Fluid-Base Interface. Mathematical Methods in the Applied Sciences, 30, 2037-2053.
https://doi.org/10.1002/mma.913

[16]   Hill, J.M. (1989) Similarity Solution for Nonlinear Diffusion—A New Integration Procedure. Journal of Engineering Mathematics, 23, 141-155.
https://doi.org/10.1007/BF00128865

[17]   Hill, D.L. and Hill, J.M. (1990) Similarity Solutions for Nonlinear Diffusion—Further Exact Solutions. Journal of Engineering Mathematics, 24, 109-124.
https://doi.org/10.1007/BF00129869

[18]   Momoniat, E. and Abelman, S. (2004) An Investigation into the Spreading of a Thin Liquid Drop under Gravity on a Slowly Rotating Disk. International Journal of Non-Linear Mechanics, 39, 265-270.
https://doi.org/10.1016/S0020-7462(02)00173-7

[19]   Momoniat, E., Mason, D.P. and Mahomed, F.M. (2001) Non-Linear Diffusion of an Axisymmetric Liquid Drop: Group Invariant Solution and Conservation Law. International Journal of Non-Linear Mechanics, 36, 879-885.
https://doi.org/10.1016/S0020-7462(00)00051-2

[20]   Momoniat, E., Myers, T.G. and Abelman, S. (2005) New Solutions for Surface Tension Driven Spreading of a Thin Film. International Journal of Non-Linear Mechanics, 40, 523-529.
https://doi.org/10.1016/j.ijnonlinmec.2004.07.017

[21]   Acheson, D.J. (1990) Elementary Fluid Dynamics. Clarendon Press, Oxford, 238-251.

[22]   Gillespie, R.P. (1959) Integration. Oliver and Boyd, Edinburgh, 90-95 and 113-116.

[23]   Ibragimov, N.H. and Anderson, R.L. (1994) Apparatus of Group Analysis. In: Ibragimov, N.H., Ed., CRC Handbook of Lie Group Analysis of Differential Equations, Vol. 1: Symmetries, Exact Solutions and Conservtion Laws, CRC Press, Boca Raton, 12-14.

[24]   Barnard, S. and Child, J.M. (1936) Higher Algebra. MacMillan, London.

[25]   Nayfeh, A.H. (1981) Introduction to Perturbation Techniques. John Wiley and Sons, New York, 39-43.

 
 
Top