Scaling

To scale the variables we define a new vertical coordinate \(\eta = y/\delta\), where \(\delta\) is the approximate boundary layer thickness, so that:

\[ \eta = y \sqrt{\frac{U_\infty}{\nu x}} \]

The continuity equation can be integrated by introducting a stream function \(\psi\), which can be written as:

\[ \psi = \sqrt{\nu x U_\infty} \,\, f(\eta), \]

where \(f(\eta)\) is the dimensionless stream function.

The \(u\) velocity components can be written as

\[ u = \frac{\partial \psi}{\partial y} = \frac{\partial \psi}{\partial \eta} \frac{\partial \eta}{\partial y} = U_\infty f'(\eta)\]

where the prime represents differentiation wrt to \(\eta\). Similarly the \(v\) component is

\[ v = - \frac{\partial \psi}{\partial x} = \frac{1}{2} \sqrt{\frac{\nu U_\infty}{x}} (\eta \, f' - f)\]

Other required terms can be further derived, including \(\partial u/\partial x\), \(\partial u/\partial y\), and \(\partial^2 u/\partial y^2\):

\[ \frac{\partial u}{\partial x} = U_\infty \, f''(\eta) \left(-\frac{1}{2} \right) \frac{\eta}{x} = -\frac{U_\infty \eta}{2 x} \, f''(\eta) \]

\[ \frac{\partial u}{\partial y} = U_\infty \,f''(\eta) \frac{\partial \eta}{\partial y} = U_\infty \sqrt{\frac{U_ \infty}{\nu x}} \, f''(\eta) \]

\[ \frac{\partial^2 u}{\partial y^2} = U_\infty \sqrt{\frac{U_\infty}{\nu x}} \left( f'''(\eta) \sqrt{\frac{U_\infty}{\nu x}} \right) = \frac{U_\infty^2}{\nu x} \, f'''(\eta) \]

Substituting all terms into the \(x\) momentum equation gives

\[ -\frac{U_\infty^2}{2 x} \eta \, f' \, f'' + \frac{U_\infty^2}{2 x} (\eta \, f' - f) \, f'' = \frac{U_\infty^2}{x} \, f'''\]

which can be simplified to

\[ f \, f'' + 2\, f''' = 0 \]

or, to be consistent with the later numerical approach,

\[ f''' + \frac{1}{2} \, f \, f'' = 0 \]

Boundary conditions

The boundary conditions correspond to

\[ f(\eta=0) = 0, \,\,\, f'(\eta=0) = 0 \]

and

\[ f'(\eta = \infty) = 1 \]

The latter BC describes the fact that at large \(y\) the value must approach the free stream velocity \(U_\infty\).

The resulting differential equation is both nonlinear and of the third order, thus the 3 BCs are sufficient to determine the solution. However, there is no exact analytical solution, and must therefore be solved either approximately (e.g. by a series expansion method as done by Blasius) or numerically (as done by Howarth in 1938⁴). One possible numerical approach will be explored in the next section.

4th order Runge-Kutta scheme⁶

The classic Runge-Kutta method for numerically solving a differential equation involves creating a weighted average of slopes and using that combined with the previous value to estimate each new value (see Figure 2). When run from an initial condition, this can be used to solve the differential equation. Classically this is thought of as an integration in time (in oceanography, I hear the terms “Runge-Kutta” most often when listening to talks by modelers who use their model velocity output to simulate drifters), but it can be used for any generic DE system.

Figure 2: Runge-Kutta method, from Wikipedia

To implement the RK4 algorithm, we create a function that implements the weighted average algorithm based on the input function fz\(=f(z)\), it’s first derivative fp\(=f'(z)\), and the step size dz.

rk4 <- function(fz, fp, dz) {
    k1 <- fp(fz)
    k2 <- fp(fz+dz/2*k1)
    k3 <- fp(fz+dz/2*k2)
    k4 <- fp(fz+dz*k3)
    return(fz + dz/6*(k1+2*k2+2*k3+k4))
}

In order to apply the RK4 method to the Blasius equation, which is a third-order DE, we have to replace it with three first-order DEs, e.g.

\[ F' = G \]

\[ G' = H \]

\[ H' = F''' = -\frac{1}{2} F H \]

To do this, we create a function for the first derivatives so that we can simultaneously solve the coupled system:

fp <- function(f) {
    return(c(f[2], f[3], -0.5*f[1]*f[3]))
}

We can use the initial conditions and set up the arrays for the solution, according to \(F(0) = 0\), \(F'(0) = G(0) = 0\), and \(F'(\infty) = G(\infty) = 1\). However, we also need a condition for \(F''(0) = H(0)\), which is unknown. This is where the shooting method comes in; we make a guess for \(H(0)\), solve the system, and check whether the far-field boundary condition \(F'(\infty) = 1\) is satisfied. \(H(0)\) can then be adjusted to achieve the far-field condition.

Choosing a random value of \(H(0) = 0.5\), the solution produces:

H0 <- 0.5

dz <- 0.05
Z <- 10
N <- Z/dz
z <- seq(0, Z-dz, dz)
F <- G <- H <- rep(0, N)
H[1] <- H0
for (i in 1:(N-1)) {
    fz <- c(F[i], G[i], H[i])
    fz2 <- rk4(fz, fp, dz)
    F[i+1] <- fz2[1]
    G[i+1] <- fz2[2]
    H[i+1] <- fz2[3]
}
plot(z, G)
grid()
abline(h=1)

Remember that \(G(\infty) = 1\), so it’s clear that our guess for \(H(0)\) is too high. Trial and error suggests that an appropriate value is \(H(0) \sim 0.33\), however we can easily do better than this by wrapping the solver in a loop and iterating based on an error tolerance:

Hmin <- 0
Hmax <- 1
err <- 1
counter <- 1
dz <- 0.05
Z <- 10
N <- Z/dz
z <- seq(0, Z-dz, dz)
HH <- NULL
while (abs(err) > 1e-6) {
    H0 <- (Hmin + Hmax)/2
    HH[counter] <- H0
    F <- G <- H <- rep(0, N)
    H[1] <- H0
    for (i in 1:(N-1)) {
        fz <- c(F[i], G[i], H[i])
        fz2 <- rk4(fz, fp, dz)
        F[i+1] <- fz2[1]
        G[i+1] <- fz2[2]
        H[i+1] <- fz2[3]
    }
    err <- tail(G, 1) - 1
    if (err < 0) Hmin <- H0 else Hmax <- H0
    ## plot(z, G)
    ## abline(h=1)
    counter <- counter + 1
}

par(mfrow=c(2, 1), mar=c(3, 3, 1, 1), mgp=c(1.5, 0.5, 0)) 
plot(z, G); grid()
plot(HH, type='l', ylim=c(0, 1), xlab='Iteration', ylab='H0 guess'); grid()
mtext(paste0('H0 = ', format(tail(HH, 1), digits=6)), side=3, line=-2, adj=1, cex=2)

References