sandbox/easystab/orr_sommerfeld_squire.m

    The Orr-Sommerfeld and Squire equation

    Here we do just as in poiseuille_uvp.m, except that instead of using the primitive variables formulation (u,v,p), we use the Orr-Sommerfeld and Squire equations, with the variables v and w. We test it on the Poiseuille flow.

    clear all; clf

N=200; % number of grid points
alpha=1; % wavenumber in x
beta=0; % wavenumber in z
Re=10000; % Reynolds number

% differentiation matrices 
[y,DM] = chebdif(N,4); 
d.y=DM(:,:,1);    
d.yy=DM(:,:,2);    
d.yyyy=DM(:,:,4);    

I=eye(N);Z=zeros(N,N);II=blkdiag(I,I);
d.x=i*alpha*I;
d.z=i*beta*I;

% base flow and derivatives
u=1-y.^2;
up=-2*y;
upp=-2;

% laplacian
k2=alpha^2+beta^2;
lap=(d.yy-k2*I); 
lap2=(d.yyyy-2*k2*d.yy+k2^2*I); 

% system matrices
LOS=-diag(u)*lap*d.x+diag(upp)*d.x+lap2/Re;     
LSQ=-diag(u)*d.x+lap/Re;
LC= -diag(up)*d.z;

A=[LOS,Z;LC,LSQ]; 
E=[lap,I];

    Boundary conditions

    They are: zero at the walls for v and w (no-slip) and zero derivative at the walls for v (that comes from the no-slip condition on u). The contraint matrix is thus \displaystyle C=\left(\begin{array}{c} I|_{-1}&0\\ I|_{1}&0\\ 0&I|_{-1}\\ 0&I|_{1}\\ \partial_y|_{-1}&0\\ \partial_y|_{1}&0\\ \end{array}\right) such that the boundary conditions are implemented \displaystyle Cq=0 where q is the state vector \displaystyle q=\left(\begin{array}{c} v\\ w \end{array}\right) 

    % boundary conditions
v=1:N; w=v+N; % location vectors
loc=[v([1,2,N-1,N]),w([1,N])];
C=[II([v([1,N]),w([1,N])],:); ... % Dirichlet on v and w
    d.y([1,N],:), Z([1,N],:)];    % Neumann on v

A(loc,:)=C;
E(loc,:)=0;
    
% computing eigenmodes 
[U,S]=eig(A,E);
s=diag(S);  [t,o]=sort(-real(s)); s=s(o); U=U(:,o);
rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[];

% loading the scanned figure for comparison and plotting
a=imread('poiseuille_spectra.png'); ss=size(a);
xx=linspace(0,1,ss(2));
yy=linspace(0,-1,ss(1));
image(xx,yy,a); axis xy
hold on

plot(-imag(s)/alpha,real(s),'ro')
grid on
axis([0,1,-1,0.1]);
xlabel('wave speed'); ylabel('exponential growth rate')

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','orr_sommerfeld_squire.png');
    The spectrum

    The spectrum

    Exercices/Contributions

    • Please write down the equations and how to obtain them in the comments
    • Please test it for a case with non-zero \beta (an oblique wave)
    • Please test it for the Blasius boundary layer
    • Please do a quantitative comparison of the accuracy with the primitive variables formulation poiseuille_uvp.m.