Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
Modelling the flow of a thin film of fluid over a
horizontal substrate
Introduction
The flow of thin films of viscous fluids, such as honey or syrup, can be effectively modelled
by identifying the dominant forces driving the flow and applying the approximations of
lubrication theory, which hold as long as the thickness of the layer of fluid is much smaller
than its horizontal extent. A schematic diagram of a typical thin film of viscous fluid,
of thickness H, spreading under gravity over a horizontal substrate is shown in figure 1.
Such flows are also referred to as viscous gravity currents, and are common in the world
around us, ranging from the spread of honey over morning toast, to lava flows and the
flow of glacial ice sheets [1, 5, 6, 7, for example].
Assuming that vertical shear stress provides the dominant resistance to the flow,
and that inertial effects and the effects of surface tension are negligible, the flow can be
shown to be governed by the following momentum balance equation
0 = −∇p− ρgez + µ
∂2u
∂z2
, (1)
where the velocity field u = (u, 0, 0) is horizontal, and ez = (0, 0, 1) is the unit basis
vector in the z-direction. Here, the pressure p and velocity u may, in general, depend
on all spatial coordinates (x, y, z) as well as on time t, whereas the density ρ and dy-
namic viscosity µ are constants which describe the material properties of the fluid. The
gravitational acceleration is denoted by g. Axes are chosen so that the fluid flows along
the x-direction, z = 0 describes the horizontal substrate, and z = H describes the upper
surface of the thin film of viscous fluid as shown in figure 1. The viscous fluid is supplied
at a specified flux at x = 0, and its leading edge is given by x = xN (t). The setup is
fully two-dimensional, with no variation along the y-direction.
Equation (1) is a version of the so-called Stokes equations, used to describe the flow
of viscous fluids.
xN(t)
H(x,t)
x
z
inflow
Figure 1: Schematic diagram illustrating the side profile of a thin film of viscous fluid
spreading under gravity over a rigid horizontal substrate.
1
Problem 1: hydrostatic pressure
Solve for the pressure by considering the z-component of (1) and applying the boundary
condition that the pressure attains its atmospheric value p0 at the upper surface of the
viscous fluid, that is
p = p0 at z = H. (2)
Problem 2: equation for the horizontal velocity
By substituting your expression for the pressure into the x-component of (1) and as-
suming that the atmospheric pressure p0 is constant, show that the horizontal velocity
u satisfies
∂2u
∂z2
=
ρg
µ
∂H
∂x
. (3)
Problem 3: solving the Stokes equations
Noting that the right-hand side of (3) is independent of z, integrate both side of (3)
twice in z in order to obtain an expression for u in terms of two arbitrary constants, A
and B.
Problem 4: velocity field
Solve for A and B by applying the boundary condition that the upper surface is stress
free, that is,
∂u
∂z
= 0 at z = H, (4)
and that the fluid does not slip at the contact line with the substrate, that is,
u = 0 at z = 0. (5)
The latter is also known as the no-slip condition.
Problem 5: depth-integrated flux
The depth-integrated flux of fluid per unit width is given by
q =
∫ H
0
u dz. (6)
After substituting in your determined values of A and B into your expression for the
velocity, show that the depth integrated flux is given by
q = −CHk
∂H
∂x
. (7)
where C and k are constants that you should determine.
2
Problem 6: Similarity solution
Along with the equation (7) specifying the flux q, the deformation of the upper surface
is governed by the mass conservation equation
∂H
∂t
= −
∂q
∂x
, (8)
along with the source flux boundary condition
q = Q at x = 0, (9)
the vanishing frontal thickness condition
H = 0 at x = xN (t), (10)
and the global mass conservation condition
∫ xN
0
H dx = Qt. (11)
By rescaling the surface height H as
H(x, t) = αtaF (η) (12)
where
η =
x
βtb
and xN = ηNβt
b, (13)
for some constant ηN , find the values of a, b, α and β so that the governing equations
(7) and (8) reduce to the single ordinary differential equation
1
5
F −
4
5
η
dF
dη
=
1
3
d
dη
(
F 3
dF
dη
)
, (14)
and recast the boundary and integral conditions (9)–(11) in terms of only F and its
derivatives. No physical constants, such as µ, ρ, Q, g, and ηN
1 and no time variable t
should appear in these reduced boundary and integral conditions. If any of these appear
in your boundary and integral conditions, that means that your values of a, b, α and
β should be adjusted. A similar, though differently scaled, derivation can be found in
[4] for Newtonian viscous fluids and in [8] for non-Newtonian, power-law fluids. A more
general treatment of self-similar problems, and ways of obtaining similarity solutions,
can be found in [2].
Note: the solution F is known as a similarity solution and the variable η is known
as a similarity variable. Note also that the front x = xN corresponds to η = ηN under
the transformation (13).
1
ηN will actually appear, but only in your transformed version of (10) and (11).
3
Problem 7: Asymptotic solution near the front
In the neighbourhood of the front η = ηN (or, equivalently, x = xN ), the thickness of the
liquid layer is small, while its slope is large. Find an asymptotic solution by following
the steps outlined below.
(a) By changing variables so that
η = ηN − ǫX, (15)
where ǫ is a small parameter, and rescaling F as F = ǫ1/3G, determine an ordinary
differential equation forG in terms ofX only (no η should appear in your equation).
(b) Identify any terms that are small in comparison to the remaining terms, when
ǫ→ 0.
(c) Formulate an approximate equation for G by deleting the small terms identified
above. Your approximate equation should have no ǫ’s and no η’s in it; just ηN and
G as a function of X. If you have any ǫ’s in your equation, that means that you
either did not delete all the small terms, or you did not divide both sides of the
equation by the right power of ǫ.
(d) Show that the solution to this approximate equation is of the form
G ∼ γXr (16)
where γ and r are constants that you should determine.
(e) Recast your approximate solution for G in terms of F and η. This is known as an
asymptotic solution, and it is valid near the front η = ηN .
(f) For what value of ηN does your asymptotic solution for F satisfy the source flux
boundary condition (the version of (9) recast in similarity coordinates for Problem
5)?
(g) Plot the asymptotic solution in Mathematica, or another program, as a function
of η over the domain [0, ηN ] for the value of ηN found in part (f), above.
A more general overview of asymptotics methods, and ways of obtaining asymptotic
solutions to differential equations, can be found in [3].
References
Project A2
A numerical treatment of similarity solutions
in the real world: thin film flows
Introduction
Many physical phenomena exhibit so-called self-similar behaviour, in which the physical
system grows or decays in time while retaining its general shape and structure, albeit
in rescaled (self-similar) coordinates. For such systems, any initial irregularities (as
formulated via initial conditions) eventually evolve towards a self-similar regime, in which
these initial irregularies become ‘forgotten’ with time. This project will lead you through
an example of such a physical system by examining the fluid mechanics of a thin film
of viscous fluid spreading under gravity over a horizontal, rigid substrate. Such flows
are also referred to as viscous gravity currents, and are common in the world around
us, ranging from the spread of honey over morning toast, to lava flows and the flow of
glacial ice sheets [1, 5, 6, 7, for example].
The key feature behind such flows is that they involve very viscous fluids, such as
honey, or syrup, and as such, viscous forces determine the flow, and any inertial effects
are negligible. Under appropriate assumptions, including that the flow is long and thin,
and resisted dominantly by vertical shear stresses, and that inertial effects and the effects
of surface tension are negligible, the depth-integrated flux of viscous fluid per unit width
can be shown to be given by
q = −
ρg
3µ
H3
∂H
∂x
, (1)
where ρ, µ and H = H(x, t) are the density, dynamic viscosity, and thickness of the
viscous fluid, respectively, and g is the acceleration due to gravity. Here, the flow is fully
two-dimensional, with axes chosen so that the fluid flows along the x-direction, z = 0
describes the horizontal substrate, z = H describes the upper surface of the thin film
of viscous fluid, and there is no variation along the y-direction. The viscous fluid is
supplied at a specified flux at x = 0, and its leading edge, or frontal position, is given
by x = xN (t) as shown in figure 1.
xN(t)
H(x,t)
x
z
inflow
Figure 1: Schematic diagram illustrating the side profile of a thin film of viscous fluid
spreading under gravity over a rigid horizontal substrate.
1
Along with the equation (1) specifying the flux q, the deformation of the upper
surface is governed by the mass conservation equation
∂H
∂t
= −
∂q
∂x
, (2)
along with the source flux boundary condition
q = Q at x = 0, (3)
the vanishing frontal thickness condition
H = 0 at x = xN (t), (4)
and the global mass conservation condition
∫ xN
0
H dx = Qt. (5)
Problem 1: Similarity solution
By rescaling the surface height H as
H(x, t) = αtaF (η) (6)
where the similarity variable η is given by
η =
x
βtb
and xN = ηNβt
b, (7)
for some constant ηN , find the values of a, b, α and β so that the governing equations
(1) and (2) reduce to the single ordinary differential equation
1
5
F −
4
5
η
dF
dη
=
1
3
d
dη
(
F 3
dF
dη
)
, (8)
and recast the boundary and integral conditions (3)–(5) in terms of only F and its
derivatives. No physical constants, such as µ, ρ, Q, g, and ηN
1 and no time variable t
should appear in these reduced boundary and integral conditions. If any of these appear
in your boundary and integral conditions, that means that your values of a, b, α and
β should be adjusted. A similar, though differently scaled, derivation can be found in
[4] for Newtonian viscous fluids and in [8] for non-Newtonian, power-law fluids. A more
general treatment of self-similar problems, and ways of obtaining similarity solutions,
can be found in [2].
Note: the solution F is known as a similarity solution and the variable η is known
as a similarity variable. Note also that the front x = xN corresponds to η = ηN under
the transformation (7).
1
ηN will actually appear, but only in your transformed version of (4) and (5).
2
Problem 2: Asymptotic solution near the front
In the neighbourhood of the front η = ηN (or, equivalently, x = xN ), the thickness of the
liquid layer is small, while its slope is large. Find an asymptotic solution by following
the steps outlined below.
(a) By changing variables so that
η = ηN − ǫX, (9)
where ǫ is a small parameter, and rescaling F as F = ǫ1/3G, determine an ordinary
differential equation forG in terms ofX only (no η should appear in your equation).
(b) Identify any terms that are small in comparison to the remaining terms, when
ǫ→ 0.
(c) By deleting the small terms identified above, show that G satisfies the following
approximate equation
4
5
ηNG
′ =
1
3
(G3G′)′. (10)
(d) Show, by substitution, that this approximate equation has a solution of the form
G ∼ γXr (11)
where γ and r are constants that you should determine.
(e) Recast your approximate solution for G in terms of F and η, and call it Fasym.
You should have no ǫ’s in your formula for Fasym. This is known as an asymptotic
solution, and it is valid near the front η = ηN .
A more general overview of asymptotics methods, and ways of obtaining asymptotic
solutions to differential equations, can be found in [3].
Problem 3: Numerical solution – free boundary problem
Follow the steps outlined below to solve the ordinary differential equation (8), subject
to appropriate boundary conditions, in Mathematica.
(a) Near the front of the viscous fluid (that is, near η = ηN ), the slope of the solution
is negative and large, approaching −∞ as η → ηN . All numerical solvers, no
matter how good, run into problems when the desired solution becomes infinite at
any point in the domain. To eliminate such problems, instead of solving over the
entire interval [0, ηN ], you should instead solve over a shorter interval [0, ηN − ǫ],
where ǫ is a small number. You should use the asymptotic solution Fasym found
in Problem 2(e), as the solution for the interval [ηN − ǫ, ηN ] – the remaining part
of the interval. The asymptotic solution will also be useful in order to formulate
appropriate boundary conditions at η = ηN − ǫ, in order for the numerical solver
3
to be able to give you the correct numerical solution. To do so, substitute your
formula for the asymptotic solution Fasym into the boundary conditions formulated
below:
F (ηN − ǫ) = Fasym(ηN − ǫ) (12)
F ′(ηN − ǫ) = F
′
asym(ηN − ǫ) (13)
(b) Use the Mathematica built-in numerical differential equation solver NDSolve to
solve the ordinary differential equation (8) subject to the boundary conditions
(12) and (13) over the interval [0, ηN − ǫ], assuming that ηN = 1 and ǫ = 0.001,
for now. You will want to be able to change the parameters ηN and ǫ later on,
so you should keep these parameters general when typing in (8), (12) and (13)
into Mathematica, and only assign a specific value to them in NDSolve using a
substitution rule – see
(c) Once you obtain the solution, plot the surface height F (η) as a function of η over
the interval [0, ηN − ǫ], save it (right-click the plot in Mathematica, and select Save
Graphic As) and insert your plot into your write-up.
(d) Use your numerical solution obtained in part (c) above to evaluate the flux
qsim = −
1
3
F 3(0)F ′(0) (14)
at η = 0.
(e) According to one of your boundary conditions obtained for Problem 1, this flux,
given in (14), should be 1. As the answer you obtained in part (d), above, is not
equal to 1, this means that the parameter ηN should be adjusted. Recall we have
just arbitrarily chosen ηN = 1 in part (c), above, for the sake of obtaining just
some numerical solution. You should now manually adjust (or write a short script
to programatically adjust) the value of ηN so that the boundary condition
qsim = −
1
3
F 3(0)F ′(0) = 1 (15)
is satisfied at η = 0, to within a tolerance of 10−3. State the value of ηN that
you have found. Plot the surface height F (η) as a function of η over the interval
[0, ηN − ǫ], save it and insert your plot into your write-up.
(f) Overlay on top of your previous plot a plot of the asymptotic solution Fasym over
the whole interval [0, ηN ]. Save this plot and add it to your write-up.
(g) Comment on how well the asymptotic solution approximates the numerical solution
over the entire interval. In which part of the interval is the fit best?
4