1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Blasius Model in Matlab

  1. Mar 8, 2016 #1
    1. The problem statement, all variables and given/known data
    Program, without any built in functions (like ODE45), a solution to the Blasius Equation in Matlab that outputs boundary layer profiles for given x values, u values, etc.

    2. Relevant equations
    2f''' + f''f = 0
    fj = fj-1 + Δη/2 * (gj + gj-1)
    gj = [Δη2/4][2gj+1/Δη2 + 2gj-1/Δη2 + fj*(gj+1 - gj-1)/2Δη]
    η = y[U/νx]1/2

    3. The attempt at a solution
    So I'm not sure how to do this in Matlab. What do I put for f? I understand I can loop through once I have j-1, j, and j+1 terms, but I guess getting there and getting them into a form that I can manipulate evades me. And then afterwards what to do. How to I relate η to f in order to create boundary layer profiles? I can't use any built in functions like ODE45, but I wouldn't even know how to use that either so I guess that's ok!
  2. jcsd
  3. Mar 8, 2016 #2


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    You can't solve any differential equation numerically unless you know some initial conditions to start with.

    The Blasius equation is discussed in more detail here:


    Doing a web search will turn up other hits.

    As far as the mechanics of the solution are concerned, you are forbidden to use canned routines like ODE45, but that doesn't mean you are exempt from using the method, just that you must code it up for Matlab. You should study Runge-Kutta type methods for solving ODEs.
  4. Mar 8, 2016 #3
    The B.C.'s we have are f(0) = 0, f'(0) = 0, and f'(∞) = 1, and that's directly relatable to η (so pretty much f(η), or at least that's what it looks like).

    I guess that math sort of confuses me than. I'm looking at this without much going on upstairs. I have a whole bunch of letters, like f and g, and I know I'm solving for f and g, but I have no clue how I can solve for f and g when the definition of them include f1, g1, etc. And that link is showing that f is a function of x and y, among other things? But I don't have an f . . .
  5. Mar 8, 2016 #4


    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    If you've never solved an ODE numerically, the notation can look pretty confusing. This seems like a difficult problem to start with, being a third order equation.

    Numerical techniques are usually designed to solve first order equations, which means you'll have to convert the Blasius equation into an equivalent system of first order ODEs in order to solve it numerically. Have you been shown how to do this?
  6. Mar 8, 2016 #5
    . . . Possibly I have learned that? . . .

    So here's the equations for fj and gj we were given:
    gj = [∆η2/4][2gj+1/∆η2 + 2gj-1/∆η2 + fj([gj+1-gj-1]/2∆η)]
    fj = ∆η[fj-1/∆η + (gj+gj-1)/2] = fj-1 + ∆η(gj + gj-1)

    So to solve that, do I need to put in random initial guess to help iterate through? Because I have no gj+1 terms . . . or any j+1 terms for that matter. And nothing is really making sense when I read through it.
  7. Mar 9, 2016 #6
    I found this Matlab code:

    % Program: blasius
    % Solve the Blasius problem
    % using 2nd order Runge-Kutta
    % Mark Blyth 2007

    clear all
    close all
    clc % Clears the command window

    % Set initial condition, step length (h)
    % and termination point

    h= 0.05; tend = 5.0;


    y1(1)=0; y2(1)=0; y3(1)=alpha;

    % Advance solution up to tend



    while t(i) < tend


    % Step from t(i-1) to t(i) using
    % the provisional step tt inbetween

    % 1. Compute f(t,y)

    f10 = y2(i-1);
    f20 = y3(i-1);
    f30 = -y1(i-1)*y3(i-1);

    % 2. Set provisional values

    yt1 = y1(i-1) + h*f10;
    yt2 = y2(i-1) + h*f20;
    yt3 = y3(i-1) + h*f30;

    % 3. Compute phase speed ft at provisional point

    ft1 = yt2;
    ft2 = yt3;
    ft3 = -yt1*yt3;

    % 4. Compute final phase speed ff

    ff1 = (1/2)*f10 + (1/2)*ft1;
    ff2 = (1/2)*f20 + (1/2)*ft2;
    ff3 = (1/2)*f30 + (1/2)*ft3;

    % 5. Set new solution value at t(i)

    y1(i) = y1(i-1) + h*ff1;
    y2(i) = y2(i-1) + h*ff2;
    y3(i) = y3(i-1) + h*ff3;

    % 6. Update time

    t(i) = t(i-1) + h;


    % Plot solution for df/dx

    hold on

    How does this make the profile? I don't see anything in there that looks familiar, and I'm trying to dissect these codes to maybe get sent on the right track. And how can you use this to create velocity profiles?
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted