boundary_layer_2D.m

A code that looks very much like jet_2D.m except for the fact that the top boundary is not a slip boundary. We build the inflow profile with a Blasius boundary layer profile to see how it behaves in 2D with the Navier-Stokes equations. The Blasius profile is computed in blasiusf.m, a function version of the code blasius.m.

See jet_2D.m and then venturi.m to learn how this works.

Dependency:

clear all; clf; format compact

% parameters 
Re=2; % Inflow reynolds number
Nx=100; % number of grid nodes in z
Ny=50; %number of grid nodes in r
Lx=100; % domain length
Ly=30; % domain height

% differentiation 
[d.x,d.xx,d.wx,x]=dif1D('fd',0,Lx,Nx,5);
[d.y,d.yy,d.wy,y]=dif1D('cheb',0,Ly,Ny);
[D,l,X,Y,Z,I,NN]=dif2D(d,x,y);

% preparing boundary conditions
NN=Nx*Ny;
l.u=(1:NN)'; l.v=l.u+NN; l.p=l.v+NN;
II=speye(3*NN);
D.lap=D.yy+D.xx;

neuploc=[l.ctl;l.ctr;l.ctr-Ny];  % where to impose the neumann condition on the pressure
p0loc=2*Ny; % where to impose zero pressure
dir=[l.left;l.top;l.bot]; % where to put Dirichley on u and w

loc=[l.u(dir); l.v(dir); l.p(p0loc); ...
    l.u(l.right); ...
    l.v(l.right); ...
    l.p(neuploc)];

C=[II([l.u([l.left;l.bot]);l.v(dir);l.p(p0loc)],:); ...     % Dirichlet on u,v,and p
   D.y(l.top,:)*II(l.u,:); ... % slip at top for u
   D.x(l.right,:)*II(l.u,:); ...   % Neuman on v at outflow
   D.x(l.right,:)*II(l.v,:); ...   % Neumann on u at outflow
   D.lap(neuploc,:)/Re*II(l.v,:)-D.x(neuploc,:)*II(l.p,:)]; % neuman constraint on pressure

% Compute Blasius velocity profile
[ybla,ubla]=blasiusf(Ly,100); plot(ubla,ybla)
uprof=interp1(ybla,ubla,y);

% initial guess
V=zeros(NN,1);
U=uprof*ones(1,Nx);
P=-X/Re; P=P-P(p0loc); % pressure zero at p0loc
q0=[U(:);V(:);P(:)];

% Newton iterations
disp('Newton loop')
q=q0;
quit=0;count=0;
while ~quit     
 
    % the present solution and its derivatives
    U=q(l.u); V=q(l.v); P=q(l.p);
    Ux=D.x*U; Uy=D.y*U; 
    Vx=D.x*V; Vy=D.y*V;
    Px=D.x*P; Py=D.y*P;

    % nonlinear function
    f=[-U.*Ux-V.*Uy+D.lap*U/Re-Px; ...
       -U.*Vx-V.*Vy+D.lap*V/Re-Py; ...
       Ux+Vy];
    
    % Jacobian 
    A=[-(spd(U)*D.x+spd(Ux)+spd(V)*D.y)+D.lap/Re, -spd(Uy), -D.x; ...
       -spd(Vx),-(spd(U)*D.x+spd(V)*D.y+spd(Vy))+D.lap/Re, -D.y; ...
       D.x, D.y, Z];
     
    % Boundary conditions 
    f(loc)=C*(q-q0);
    A(loc,:)=C;

    % plotting
    subplot 311;
    surf(X,Y,reshape(P-1,Ny,Nx)); view(2); shading interp; hold on
    
    sely=1:Ny; selx=1:6:Nx;
    ww=reshape(U,Ny,Nx); vv=reshape(V,Ny,Nx); 
    quiver(X(sely,selx),Y(sely,selx),ww(sely,selx),vv(sely,selx),'k');
    axis([0,Lx,0,Ly]);
    xlabel('z'); ylabel('r'); title('Pressure P'); grid off;hold off
    
    subplot 312;
    surf(X,Y,reshape(U,Ny,Nx)); view(2); shading interp; 
    xlabel('x'); ylabel('y'); title('horizontal velocity U'); grid off
    
    subplot 313;
    surf(X,Y,reshape(V,Ny,Nx)); view(2); shading interp; 
    xlabel('x'); ylabel('y'); title('vertical velocity V'); grid off
    drawnow
    
    % convergence test  
    res=norm(f);
    disp([num2str(count) '  ' num2str(res)]);
    if count>50|res>1e5; disp('no convergence');break; end
    if res<1e-5; quit=1; disp('converged'); continue; end
    
    % Newton step
    q=q-A\f;   
    count=count+1;
end


set(gcf,'paperpositionmode','auto')
print('-dpng','-r80','boundary_layer_2D.png')

The screen output for the convergence:

Blasius: converged
Newton loop
0  33.8846
1  0.067715
2  0.0011648
3  2.9002e-07
converged

And the figure, where we see how the boundary layer diffuses in the vertical direction

The figure