Boundary layers are present any time a viscous fluid (such as air) is in contact with a surface and there is relative motion between them. Perhaps the simplest example of a boundary layer is flow over a flat-plate at 0 angle of attack. For flat plate flow, you can start with the full Navier-Stokes equations and simplify them to a 3rd-order nonlinear ODE. This is often referred to as the Blasius equation (and its solution as the Blasius boundary layer) after the guy who first discovered it:
$$f^{\prime\prime\prime} + \frac{1}{2}ff^{\prime\prime} = 0.$$
The boundary conditions are no penetration:
$$f(0) = 0,$$
no slip:
$$f^\prime(0) = 0,$$
and the far field must match the freestream flow:
$$\lim_{\eta\to\infty} \:f^\prime(\eta) = 1.$$

The purpose of this post is not to detail the derivation but to show how you can obtain solutions to this ODE using a shooting method and Matlab.

It is the last boundary condition that makes solving this equation a little difficult. We would like to start at $$\eta=0$$ and integrate by marching towards the freestream. However to do this we must chose a value of $$f^{\prime\prime}(0)$$. Unless we choose this value correctly the far-field boundary condition will not be satisfied.

One solution to this is to make a guess for $$f^{\prime\prime}(0)$$, solve the ODE, check the boundary condition, and repeat using the bisection method until the boundary condition is satisfied to within a given tolerance. This is the method implemented in the script below:

function [eta, f] = blasius(eta_max)
% Initial guesses for f''(0), these values
% need to bracket the true value
a = 0.3;
b = 1.0;

% The boundary layer only extends to around
% eta = 5, we don't want to integrate more
% than necessary (also this sometimes
% becomes unstable if this value is set too high).
% Also note that this is the height where the
% freestream boundary condition is enforced.
eta_max_ode = 20;

% Maximum iterations
max_iter = 100;

% Iteration counter
count = 0;

% Tolerance for boundary condition
epsilon = 1e-10;

% Value greater than epsilon to start it off
error = 1e15;
while(abs(error) > epsilon && count < max_iter)
% Solve the ODE using ode45
sigma = (a+b)/2;
f0 = [0, 0, sigma];
[eta_ode, f_ode] = ode45(@BlasiusFunc, [0, eta_max_ode], f0);

% Perform bisection method
error = f_ode(end,2) - 1;
if(error > 0)
b = sigma;
else
a = sigma;
end

% Increment the counter
count = count + 1;
end

f = f_ode;
eta = eta_ode;

% Adjust the ends of f and eta to respect eta_max
if(eta_max > eta_max_ode)
f = [f; [(eta_max - eta_max_ode + f(end,1)), 1, 0]];
eta = [eta; eta_max];
elseif(eta_max < eta_max_ode)
f = f(eta < eta_max,:);
eta = eta(eta < eta_max);

eta = [eta; eta_max];
f_interp = interp1(eta_ode, f_ode, eta_max);
f = [f; f_interp];
end
end

function [df] = BlasiusFunc(eta, f)
df = zeros(3,1);
df(1) = f(2);
df(2) = f(3);
df(3) = -f(1)*df(2)/2;
end


And here are some plots of the solutions:
$$f$$:

$$f^{\prime}$$:

$$f^{\prime\prime}$$:

1. Ferdousi Ema 2017-08-11

Hello brother,
I used to run the code that u gave above
but it showed a error
“” ??? Input argument “eta_max” is undefined.

Error in ==> blasius2 at 48
if(eta_max > eta_max_ode)
“”
now what can i do ?

• jason 2018-09-25

The code above is just a function. To use it save it to a file called blasius.m and call it from either the command window or another script as:
[eta, f] = blasius(10);

You can then plot it with:
figure; plot(f(:,2), eta); xlabel('$$\frac{u}{U_\infty}$$', 'interpreter', 'latex'); ylabel('$$\eta$$', 'interpreter', 'latex');

2. SHIBU CLEMENT 2017-10-24

Thanks Jason for posting the Matlab code for solving Blasius equation using Runge-Kutta and Bisection methods. Instead of bisection method if we use Newton-raphson method, what changes we have to make it in the code?

Thanks

• jason 2018-09-25

Interesting idea! The biggest challenge will be to find the derivative of R = f'(\eta_{max}) - 1 with respect to f''(0). Once you find that you can apply simple 1D Newton method to find the f''(0) that makes R zero.

3. Golam Mortuja 2018-08-14

the program is not running in MatLab 2017a. Is the erro=1e15 correct(in line 25)?

• jason 2018-09-25

Hi Golam, I just ran the code in MATLAB 2017b without errors. What problem are you seeing?