Numerical solution of PDE:s, Part 10: The thin-film equation

Earlier, I showed how to solve the 1D and 2D versions of the complex Ginzburg-Landau equation, which is an example of a nonlinear partial differential equation, and which had to be linearized for solution with implicit differencing, meaning that the matrix in the linear system was different on each timestep.

Another nonlinear PDE is the so-called thin film equation, which in 2D form reads

2dthinfilm-small.gif

Here the function h(x,y,t) describes the local thickness of a film of viscous liquid located on top of a solid surface described by the xy-plane. The parameter \gamma is the surface tension of the liquid-gas interface and \mu is the viscosity of the liquid.

An unusual thing about this equation is that it’s fourth order in the spatial coordinates, while most equations in physics are second order DE:s. Some equations of continuum mechanics describing elastic bending of cylinders and plates are also fourth order, so this is not the only example.

In many cases, a corresponding equation with only one spatial coordinate is enough for describing thin film physics, and then the graph of the solution can be though of as depicting an intersection of the film at a single value of y-coordinate.

1dthinfilm-small

When discretizing this equation, we must note that the factor h^3 has to be treated explicitly to get a linear system of equations, just like what had to be done with the |A|^2 in the CGLE. Also, we will set \gamma/(3\mu ) = 1 to make the equation dimensionless. One correct way to discretize this equation leads to the system

thinfilm-discrete.gif

where the two index object \alpha_{j}^{i} is

alpha-def-small

Note that now the linear system is not tridiagonal, but heptadiagonal due to the higher order derivatives. Solution of the equation with this differencing scheme and a Gaussian initial condition h(x,0) is done with the following R language code.

library(graphics) #load the graphics library needed for plotting

lx <- 10 #length of the computational domain
lt <- 10 #length of the simulation time interval
nx <- 150 #number of discrete lattice points
nt <- 150 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- lt/nt #length of timestep

psi = c(1:nx) #array for the function h values
sol = c(1:nx)

kappa = dt/(4*dx*dx*dx*dx)

for(j in c(1:nx)) {
psi[j] = exp(-(j*dx-5)*(j*dx-5))
sol[j] = psi[j]
}

xaxis <- c(1:nx)*dx #the x values corresponding to the discrete lattice points
alpha <- c(1:nx)

A = matrix(nrow=nx,ncol=nx) #matrix for forward time evolution

for (m in c(1:nt)) { #main time stepping loop

for(j in c(1:nx)) {
alpha[j] = kappa*psi[j]*psi[j]*psi[j]
}

for(j in c(1:nx)) {
for(k in c(1:nx)) {
A[j,k]=0
if(j==k) {
if(j!=nx && j!=1) {
A[j,k] = 1 + 2*alpha[j+1] + 2*alpha[j-1] #diagonal elements
}
if(j==1) {
A[j,k] = 1 + 2*alpha[j+1]
}
if(j==nx) {
A[j,k] = 1 + 2*alpha[j-1]
}
}

if(j==k+1 && j!=1) {
A[j,k] = -alpha[j-1] #off-diagonal elements
}

if(j==k-1 && j!=nx) {
A[j,k] = -alpha[j+1]
}

if(j==k+2 && j<nx) {
A[j,k] = -2*alpha[j+1]
}

if(j==k-2 && j>1) {
A[j,k] = -2*alpha[j-1]
}

if(j==k+3 && j<nx) {
A[j,k] = alpha[j+1]
}

if(j==k-3 && j>1) {
A[j,k] = alpha[j-1]
}

}
}

for(l in c(1:nx)) {
psi[l] = sol[l]
}
sol <- solve(A,psi) #solve the system of equations

jpeg(file = paste("plot_",m,".jpg",sep=""))
plot(xaxis,sol,xlab = "position (x)",ylab="h(x,t)",ylim=c(0,1),pch='.')
title(paste("h(x,t) at t = ",round(m*dt,digits=2)))
lines(xaxis,sol)
dev.off()

}

dev.off()

An animation of the solution looks like this.

The solutions of this equation have the property that the graph settles into the shape of a downward opening parabola when time proceeds. A problem with this is that the contact angle, in which the liquid surface approaches the solid surface (x-axis) at the final equilibrium, can be anything between 0 and 90 degrees depending on the relative width and height of the initial Gaussian. In a real liquid-solid system, the final contact angle depends on the surface tensions of both the liquid-solid and the liquid-gas interfaces, as described by the Young equation.

To create an equation that can model equilibrium contact angles appropriately, we add a disjoining pressure term \Pi (h) in the equation, as here:

thinfilm-disjoin-small

One form of the \Pi -term that works is

disjoin-def-small

where the h_* is a precursor film thickness.

contact-angle.jpg

The idea behind the precursor film is that even when we have a liquid drop or puddle surrounded by apparently dry solid surface, there is actually a very thin adsorbed layer of liquid molecules of that dry area (the molecules get there by evaporating from the liquid surface and reattaching on the solid). So, in a simulation where we include the disjoining pressure, we need to use an initial condition that is a Gaussian with an added constant equal to the precursor thickness:

h(x,0) = \exp \left[-b(x-x_0 )^2 \right] + h_*

In here we will set the values n = 5 and m = 2 in the disjoining pressure term, and set the precursor film thickness into the value 0.01. The disjoining pressure term can be treated explicitly at the same time as we use implicit differencing for the rest of the equation – we subtract, on the RHS of the discretized equation a term D_{j}^i ,

discrete-disjoin

defined by

dij

where

cij

and \Pi_{j}^i is the disjoining pressure evaluated at the discrete points. Note the use of a central finite difference instead of a one-sided difference when calculating the derivatives – doing otherwise is likely to make the simulation crash. The term D_{j}^{i} only affects the right-hand side vector of the linear system \mathbf{Ax} = \mathbf{b} that we solve on each timestep. A code that solves the new equation for the pre-factor value B=0.1 is shown next.

library(graphics) #load the graphics library needed for plotting

lx <- 10 #length of the computational domain
lt <- 5 #length of the simulation time interval
nx <- 150 #number of discrete lattice points
nt <- 5000 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- lt/nt #length of timestep

prec = 0.01 #precursor film thickness

psi = c(1:nx) #array for the function h values
sol = c(1:nx)

disjoin1 = c(1:nx)
disjoin2 = c(1:nx)
disjoin3 = c(1:nx)

kappa = dt/(4*dx*dx*dx*dx)

for(j in c(1:nx)) {
psi[j] = exp(-(j*dx-5)*(j*dx-5))+prec
sol[j] = psi[j]
}

xaxis <- c(1:nx)*dx #the x values corresponding to the discrete lattice points
alpha <- c(1:nx)

A = matrix(nrow=nx,ncol=nx) #matrix for forward time evolution

for (m in c(1:nt)) { #main time stepping loop

disjoin1[1] = 0
disjoin1[nx] = 0
for(j in c(1:nx)) {
disjoin1[j] = 0.1*((prec/psi[j])^5 - (prec/psi[j])^2)
}

for(j in c(2:(nx-2))) {
disjoin2[j] = psi[j]*psi[j]*psi[j]*(disjoin1[j+1] - disjoin1[j-1])/(2*dx)
}

for(j in c(2:(nx-2))) {
disjoin3[j] = (disjoin2[j+1]-disjoin2[j-1])/(2*dx)
}

disjoin3[1] = 0
disjoin3[2] = 0
disjoin3[nx-1] = 0
disjoin3[nx] = 0

for(j in c(1:nx)) {
alpha[j] = kappa*psi[j]*psi[j]*psi[j]
}

for(j in c(1:nx)) {
for(k in c(1:nx)) {
A[j,k]=0
if(j==k) {
if(j!=nx && j!=1) {
A[j,k] = 1 + 2*alpha[j+1] + 2*alpha[j-1] #diagonal elements
}
if(j==1) {
A[j,k] = 1 + 2*alpha[j+1]
}
if(j==nx) {
A[j,k] = 1 + 2*alpha[j-1]
}
}

if(j==k+1 && j!=1) {
A[j,k] = -alpha[j-1] #off-diagonal elements
}

if(j==k-1 && j!=nx) {
A[j,k] = -alpha[j+1]
}

if(j==k+2 && j<nx) {
A[j,k] = -2*alpha[j+1]
}

if(j==k-2 && j>1) {
A[j,k] = -2*alpha[j-1]
}

if(j==k+3 && j<nx) {
A[j,k] = alpha[j+1]
}

if(j==k-3 && j>1) {
A[j,k] = alpha[j-1]
}

}
}

for(l in c(1:nx)) {
psi[l] = sol[l]
}
sol <- solve(A,psi-disjoin3) #solve the system of equations
for(l in c((nx-20):nx)) { # remove the risk of boundary effects on the right end of the domain
psi[l]=prec
sol[l]=prec
}

if(m%%10 == 0) {
jpeg(file = paste("plot_",m,".jpg",sep=""))
plot(xaxis,sol,xlab = "position (x)",ylab="h(x,t)",ylim=c(0,1),pch='.')
title(paste("h(x,t) at t = ",round(m*dt,digits=2)))
lines(xaxis,sol)
dev.off()
}

}

In the next animation, the solution curves for pre-factor values B=0.1 (red curve) and B=0.01 (black curve) are shown in the same graph.

In the animation, it is apparent that a larger value of B leads to a larger contact angle at equilibrium. Actually, it can be shown that the contact angle \theta depends on the values of the parameters as

disjoin-contact-small.gif

More information about the thin film equation and its solutions can be found on the Wiki page here, and on the NJIT department of applied mathematics homepage. When the effects of gravitation or surface tension gradients are added in the TFE, many kinds of interesting pattern formation effects can happen just like in our previous example of a nonlinear PDE, the Ginzburg-Landau equation. If the liquid film consists of a mixture of many, possibly volatile liquids, complicated multiphysics problems involving fluid dynamics, evaporation, heat transfer and chemical kinetics, all at the same time, are obtained.

 

Advertisements

Nondimensionalization and characteristic scales in physics

1 Introduction

When doing calculations in computational physics, it is usually necessary to convert equations to a nondimensional form, as a computer almost always works with dimensionless numbers. The most simple way to do this would probably be to equate the SI units (kilogram, meter, second…) to 1. This choice of scale, however, is not always appropriate as the interesting features of the dynamics of the system can take place at lengths and time intervals that are very small or large compared to these units. Plots and tables of the computational results are easier to read when the values in them are numbers that don’t differ from 1 by too many orders of magnitude. Examples of different choices of scales are the use of natural units (where we set c = 1 and ћ=1) in high-energy physics, and atomic units in quantum chemistry, where an appropriate choice of length unit is the Bohr and that of energy is Hartree.

2 Simple harmonic motion

As an example of a classical mechanical system with a characteristic timescale, consider a 1D Hookean oscillator with spring constant k and mass M. The equation of motion for the single position variable x(t) is

 M \frac{d^2 x}{dt^2} = -kx .           (1)

The obvious choice of typical timescale in this case is of course the period of oscillation

T = 2 \pi \sqrt{\frac{M}{k}} ,           (2)

and the corresponding nondimensional time variable would be \tilde{t}=t/T. This scale could also be found by simple dimensional analysis, forming a product of powers of the parameters M and k, and requiring that the product has dimensions of time:

[M]^{\alpha}[k]^{\beta} = kg^\alpha \times \frac{kg^\beta}{s^{2\beta}}=kg^{\alpha + \beta}s^{-2\beta} = s         (second) .         (3)

Solving this for alpha and beta gives us a characteristic time M^{1/2}k^{-1/2}, which is the period T divided by 2π and is of the same order of magnitude as T. If we want to integrate the equation of motion numerically with the finite-difference method, we would probably choose the discrete time step to be at least about 10-15 times smaller than T, and when plotting the trajectory x(t) we would usually choose a range of t-coordinate that is only a few periods long, no matter what the initial position x(0) and velocity x’(0) of the oscillating body are.

On the other hand, the motion of the oscillator described above does not have a general characteristic length scale. If we have two trajectories where in the first one the total energy (kinetic + potential) of the oscillator is 1 Joule and in the second one it’s 100 Joules, it is not appropriate to use the same range of x-coordinate when plotting the functions x(t), because in the latter case the amplitude of the trajectory is 10 times larger than in the former. This fact can also be seen from a simple dimensional analysis, by noting that there is no product of powers of M and k that has dimension of length:

[M]^{\alpha}[k]^{\beta} \neq m  (meter), for any \alpha, \beta .

3 Linear motion with damping

As another simple mechanical system that has a characteristic timescale, consider linear 1D motion with damping (air resistance or friction). The equation of motion is

M\frac{d^2 x}{dt^2}=-b\frac{dx}{dt},              (4)

where the damping parameter b has dimensions of mass/time. Immediately we can see that the variable M/b has dimensions of time, but we don’t know its physical significance yet. Solving the DE for x(t) and its time derivative, we obtain:

x(t) = x_0 + \frac{Mv_0}{b} - \frac{Mv_0}{b}e^{-b(t-t_0 )/M}        (5)

x'(t) = v_0 e^{-b(t-t_0 )/M} = v_0 2^{-\frac{b(t-t_0 )}{M \log 2}}=v_0 2^{-\frac{t-t_0 }{t_{1/2}}}  .      (6)

where x_0 and v_0 are the position and velocity at time t = t_0, and

t_{1/2}=\frac{M\log 2}{b}     (7)

is the half-life of the velocity.

Usually when some variable decays exponentially with some half-life, be it the number of radioactive atoms in a sample or the amount of a foreign substance (medication, environmental toxin) in human body, we can say that for our purposes the variable has practically become zero after about 6-10 half-lifes. Therefore the parameter t_{1/2}, or the related parameter M/b (which we obtained by dimensional analysis) can be seen as the characteristic timescale for the motion in this system. Again, there is no single characteristic length for the system, as the distance the object travels from the initial position x_0 in some number of half-lives is linearly proportional to the initial velocity v_0.

4  A particle inside rigid box

What about mechanical systems where the motion happens in some characteristic length scale? The simplest example of this kind of behavior is 1D motion of an object that has been confined in some interval [0, L] of the x-axis by impenetrable walls at points x = 0 and x = L. This can be described by saying that the potential energy of the object is zero when 0 < x < L and it is infinite when x < 0 or x > L. In quantum mechanics this model system is called ”particle in a box”, but here we consider only the classical mechanical equivalent.

If the collisions of the object with the walls are elastic, the trajectory x(t) is a sawtooth function like the one plotted below:

sawtooth

Figure 1. The periodic classical trajectory of a 1D particle in a rigid box.

Now the interval the object moves in is always [0,L] which doesn’t depend on its initial position and velocity (as long as the initial position is inside the box and the initial velocity is non-zero). Therefore the obvious choice of length scale for the system is L. However, this system does not have a characteristic time scale, as the period of the motion is 2L/v_0 which depends on the initial conditions.

5 A system with both a characteristic length and time

How, then, would one construct a system where there exist characteristic scales in both length and time? By dimensional analysis, we can expect that the set of parameters a,b,c,... in the equation of motion should be such that some product of their powers, a^{\alpha} b^{\beta} c^{\gamma} \dots , has dimensions of time and some other similar product has dimensions of length.

One way to make such a system is to form an equation of motion that contains both a Hookean returning force and a term that depends on both position and velocity:

M\frac{d^2 x}{dt^2} = -b x \left( \frac{dx}{dt} \right)^2 - kx        (8)

Now we can deduce that the parameter b needs to have dimensions [mass][length]^{-2}. A parameter formed from M, b and k that has dimensions of length is

L = \sqrt{\frac{M}{b}}            (9)

and a characteristic time is the familiar

T = \sqrt{\frac{M}{k}},            (10)

which already appeared in the case of simple harmonic motion. From these we can form dimensionless length and time variables

\tilde{x} = \frac{x}{L}, \tilde{t} = \frac{t}{T},            (11)

and written using these variables, the equation of motion becomes

\frac{d^2 \tilde{x}}{d\tilde{t}^2} = -A x \left( \frac{d\tilde{x}}{d\tilde{t}} \right)^2 - B\tilde{x} ,          (12)

where A and B are dimensionless numbers defined by:

A = \frac{bL^2}{M}, B = \frac{k T^2}{M}.

So how do the trajectories x(t) of this system look like? The equation of motion can’t be solved analytically, but we can do a numerical solution with Mathematica, Matlab or some other program capable of that kind of calculations. Setting the values of the parameters to A=0.5 and B=1, and using the initial position \tilde{x}(0) = 0 and a set of different initial velocities x'(0) = 0, 10, 100 or 10000, the following graphs are obtained:

v1v10

v100v10000

Figure 2. Trajectories of a particle in a velocity-dependent force field for initial dimensionless velocities of 1, 10, 100 or 10000.

From the plots in Fig. 2., it is obvious that the amplitude and period of the trajectory stay in the same orders of magnitude for a very large range of initial velocities/energies of the moving object. A 10^4 – fold increase in the dimensionless velocity causes only an about 6-fold increase in amplitude and an about 4-fold increase in frequency. The amplitude seems to increase only logarithmically with increasing velocity, i.e. slower than its any fractional power.

Note, however, that if the initial kinetic energy of the particle at t = 0 is small enough, the velocity-dependent term in the equation of motion can be neglected, and the system behaves like a simple harmonic oscillator for which the amplitude is directly proportional to the energy. To see the more interesting properties of this system, as in the graphs above, the energy needs to be large enough. In physics there are many situations where a system behaves differently in different energy scales, one example being a collision of atomic nuclei which can be described with simple Coulombic repulsion at low collision energies but will involve the possibility of nuclear fusion at larger energies where the particles can overcome the repulsive electrostatic force. The classical mechanical system here is a very simple toy model of this kind of behavior.

6 Applications

Many problems of applied physics where characteristic lengths and times and dimensionless quantities appear are related to the physics of fluid flow. A familiar dimensionless number that describes fluid motion is the Reynolds number (Re), which is related to the probability of turbulence occurring in a fluid system. Others include the capillary number (denoted Ca) and Eötvös number (denoted Eo), which describe the significance of surface tension in relation to other forces.

When a way to calculate the relevant length scales for a particular system are known, it is easier to choose a sufficiently fine discretization and a sufficiently large computational domain (intervals of the x- and t-axes) when integrating equations of motion numerically.