MATHEMATICA

TUTORIAL

under the terms of the GNU General Public License (GPL)

for the Second Course. Part VI: Blasius equation

 

---

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra

Glossary

Blasius equation

---

A Blasius boundary layer (named after the German fluid dynamics physicist Paul Richard Heinrich Blasius, 1883--1970) describes the steady two-dimensional laminar boundary layer that forms on a semi-infinite plate which is held parallel to a constant unidirectional flow. Upon introducing a normalized stream function f, the Blasius equation becomes The boundary conditions are the no-slip condition: This is a nonlinear, boundary value problem (BVP) for the third order differential equation. The main obstacle in solving this problem numerically is request to evaluate f'(∞). Several strategies have been proposed in order to bypass this evaluation. One of them is to convert the given problem into the initial value problem, which was first achieved by K. Töpfer in 1912 by constructing a transformation of variables that reduces the BVP into an IVP. Actually, this leads to determination of the second derivative at the origin, f''(0), called the Blasius constant. Its numerical value was obtained by many researchers starting with K. Töpfer; however, the rigorous derivation of the Blasius constant is due to V.P. Varin in 2013 who obtained it in the form of a convergent series of rational numbers. Therefore, we assume that its value a = f′′(0) ≈ 0.33205733621 to be known (which is slightly less than 1/3). The knowledge of the Blasius constant is also necessary to evaluate the shear stress at the plate. Once this value is determined, we obtain the initial value problem, which can be solved numerically. So we solve the relevant initial value problem using different methods.

First, we start with Picard's iteration, and to achieve this we have two options: either to apply Picard's iteration directly to the third order Blasius equation directly or to convert it to a system of first-order differential equations. So we demonstrate both options, beginning with the latter.

blasius = NDSolve[{f'''[x] + 1/2*f[x]*f''[x] == 0, f[0] == 0, f'[0] == 0, f''[0] == 1/3}, f, {x, 0, 5}] g = Plot[Evaluate[{f[x], f'[x]} /. blasius], {x, 0, 5}, PlotStyle -> {Thick, {Blue, Magenta}}] txt1 = Graphics[ Text[Style["f(x)", FontSize -> 16, Black], {3.5, 2.5}]] txt2 = Graphics[ Text[Style["f'(x)", FontSize -> 16, Black], {4.3, 0.6}]] Show[txt1, txt2, g]

 

Classical Picard's iteration

In order to apply the classical Picard's iteration, we have to convert this to a system of first-order differential equations. Let

This leads to: It is not possible to specify a boundary condition at ∞ numerically, so we will have to use a large number, and verify it is "large enough". From the solution, we evaluate the derivatives at y = 0, and we have f''(0) = f₃(0).

We have to provide initial guesses for f₁, f₂ and f₃. This is the hardest part about this problem. We know that f₁ starts at zero, and is flat there (f'(0)=0), but at large y, it has a constant slope of one. We will guess a simple line of slope = 1 for f₁. That is correct at large y, and is zero at y = 0. If the slope of the function is constant at large y, then the values of higher derivatives must tend to zero. We choose an exponential decay as a guess.

Finally, we let a solver iteratively find a solution for us, and find the answer we want. The solver is in the pycse module.

import numpy as np
from pycse import bvp

def odefun(F, x):
    f1, f2, f3 = F.T
    return np.column_stack([f2,
                            f3,
                            -0.5 * f1 * f3])

def bcfun(Y):
    fa, fb = Y[0, :], Y[-1, :]
    return [fa[0],        # f1(0) =  0
            fa[1],        # f2(0) = 0
            1.0 - fb[1]]  # f2(inf) = 1

eta = np.linspace(0, 6, 100)
f1init = eta
f2init = np.exp(-eta)
f3init = np.exp(-eta)

Finit = np.column_stack([f1init, f2init, f3init])

sol = bvp(odefun, bcfun, eta, Finit)
f1, f2, f3 = sol.T

print("f''(0) = f_3(0) = {0}".format(f3[0]))

import matplotlib.pyplot as plt
plt.plot(eta, f1)
plt.xlabel('$\eta$')
plt.ylabel('$f(\eta)$')
plt.savefig('images/blasius.png')

We get the initial second derivative to be f''(0) = 0.3324911095517483. For the sake of simplicity, we approximate the Blasius constant as 1/3.

Now we use Picard iteration procedure starting with the initial conditions. Since the derivative operator is unbounded, we apply its inverse to reduce the given initial value problem into fixed point form:

From this fixed point equation, we derive the iteration procedure The initial approximation is just the initial conditions: The first iteration becomes The second one is The third one is For next one, we have So ur fifth approximation becomes

 

Picard's iteration applied to the third order equation

We use the fixed point formulation of the given initial value problem:

This leads to the following iteration procedure: with the following initial values: For simplicity, we employ the Blasius constant to be 1/3. Next iteration results We proceed with Mathematica. Now we plot the fifth approximation along with the numerical solution obtained by Mathematica.

Plot[{Evaluate[{f[x]} /. blasius], y5[x], f5[x]}, {x, 0, 2.4}, PlotStyle -> {Thick, {Blue, Orange, Magenta}}]

Therefore, we see that the two Picard's iteration procedures do not coincide and provide two different approximations. Moreover, the Picard's iteration applied directly to the third order equation converges much faster than the iteration based on reduction to the first order system of equations.

 

Power Series Method

 

Adomian Decomposition Method

text goes here

 

References

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions