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}\):**

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 ?

please tell me about the problem

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');

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

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.the program is not running in MatLab 2017a. Is the erro=1e15 correct(in line 25)?

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