blasiusf.m

Simply the code for the Blasius boundary layer profile blasius.m, made into a function.

Input:

L is the domain height (typically 10)
N is the numbre of grid points (Typically 100)

Output

y are the grid poins
u the velocity profile on this grid

function [y,u]=blasiusf(L,N);

% differentiation matrices
scale=-2/L;
[y,DM] = chebdif(N,3); 
dy=DM(:,:,1)*scale;    
dyy=DM(:,:,2)*scale^2;    
dyyy=DM(:,:,3)*scale^3;
y=(y-1)/scale;  	
I=eye(N); Z=zeros(N,N);

% initial guess from order 4 Polhausen profile
scalepol=6*sqrt(35/37);
eta=y/scalepol; deta=dy*scalepol; % rescaled vertical coordinate
u0=1-(1-eta).^3.*(1+eta);
u0(eta>1)=1;

A=deta; A(1,:)=I(1,:); u0(1)=0; % set up integration 
g0=A\u0; % compute integral

sol=g0;

% Newton iterations
quit=0;count=0;
while ~quit
    
    % the present solution and its derivatives
    g=sol; gy=dy*g; gyy=dyy*g; gyyy=dyyy*g;
    
    % nonlinear function
    f=2*gyyy+g.*gyy;
    
    % analytical jacobian
    A=2*dyyy+diag(g)*dyy+diag(gyy)*I;

    % Boundary conditions
    loc=[1,2,N];
    C=[I(1,:); dy(1,:); dy(N,:)];
    f(loc)=C*g-[0;0;1];
    A([1,2,N],:)=C;
    
    % convergence test
    res=norm(f);
    disp([num2str(count) '  ' num2str(res)]);
    if count>50|res>1e5; disp('Blasius: no convergence'); break; end
    if res<1e-5; quit=1; disp('Blasius: converged'); continue; end
    
    % Newton step
    sol=sol-A\f;
    count=count+1;
end

% the velocity profile
u=dy*g;