Blasius Model Solution in Matlab | Boundary Layer Profiles without ODE45

[Bluestribute]

Homework Statement

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.

Homework Equations

2f^''' + f^''f = 0
f_j = f_j-1 + Δη/2 * (g_j + g_j-1)
g_j = [Δη²/4][2g_j+1/Δη² + 2g_j-1/Δη² + f_j*(g_j+1 - g_j-1)/2Δη]
η = y[U/νx]^1/2

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!

 

[SteamKing]

Homework Statement

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.

Homework Equations

2f^''' + f^''f = 0
f_j = f_j-1 + Δη/2 * (g_j + g_j-1)
g_j = [Δη²/4][2g_j+1/Δη² + 2g_j-1/Δη² + f_j*(g_j+1 - g_j-1)/2Δη]
η = y[U/νx]^1/2

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!

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:

https://en.wikipedia.org/wiki/Blasius_boundary_layer

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.

 

[Bluestribute]

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 . . .

 

[SteamKing]

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 . . .

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?

 

[Bluestribute]

. . . Possibly I have learned that? . . .

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

So to solve that, do I need to put in random initial guess to help iterate through? Because I have no g_j+1 terms . . . or any _j+1 terms for that matter. And nothing is really making sense when I read through it.

 

[Bluestribute]

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;

alpha=0.4696;

y1(1)=0; y2(1)=0; y3(1)=alpha; %%%%%%%%%
% Advance solution up to tend
%%%%%%%%%

t(1)=0;

i=1;

while t(i) < tend

i=i+1;

%%%
% 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;

end%%
% Plot solution for df/dx
%%

plot(t,y2)
hold on
plot(t,y3,'--')

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?