TS_PlanePoiseuille.m

Stability of the plane Poiseuille flow in primitive variables.

This document belongs to the easystab project, check the main page to understand the general philosophy of the project.

The present program was specifically designed as a support for the course “Introductions to hydrodynamical instabilities” for the “DET” Master’s cursus in Toulouse. For the underlying theory please see Lecture notes for chapter 7

To use this program: click “raw page source” in the left column, copy-paste in the Matlab/Octave editor, and run.

Stability of the plane Poiseuille flow in primitive variables.

This code is adapted from sandbox/easystab/poiseuille_uvp.m.

The numerical treatment of the problem is the same as in this initial program, but here the eigenvalue calculation is done inside a function “EV_Poiseuille” to allow to perform loop over the parameters.

Definition of geometry and differential operators

clear all; clf;
global y dy dyy Z I INT n
L=2;        % the height of the domain, from -1 to 1
n=100;      % the number of grid points

% differentiation and integration operators
scale=-2/L;
[y,DM] = chebdif(n,2); 
dy=DM(:,:,1)*scale;    
dyy=DM(:,:,2)*scale^2;    
y=y/scale; 
Z=zeros(n,n); I=eye(n); 
INT=([diff(y)',0]+[0,diff(y)'])/2;

Validation

The valudation is done in the program sandbox/easystab/poiseuille_uvp.m by comparing with results from Drazin & Reid and is not reproduced here

Curves lambda(k) and c(k) for various values of Re

for Re = [5772 6000 6500 7000 8000 10000 20000 50000 100000] 
 alphatab = [0.5:0.05:1.5];
 stab = [];
  for alpha = alphatab
  s = EV_Poiseuille(Re,alpha); 
  stab = [stab s];
end
    
  figure(1);
   subplot(2,1,1);
   plot(alphatab,real(stab));hold on;
   grid on;
   subplot(2,1,2);
   plot(alphatab,-imag(stab)./alphatab); hold on;
   pause(0.1); 
    
end
    
   figure(1);subplot(2,1,1);
   title('Growth rate');
   xlabel('k');ylabel('\lambda_r');    figure(1);subplot(2,1,2);
   title('Phase velocity');
   xlabel('k');ylabel('c_r');
   legend('Re = 5772', 'Re=6000', 'Re=8000','Re=10 000','Re=20 000','Re= 50 000', 'Re = 100 000');
    
   set(gcf,'paperpositionmode','auto');
   print('-dpng','-r100','PlanePoiseuille_SigmaOmegaK.png');

Growth rate and phase velocity as function of k for various Re

Construction of the Marginal stability curve

  Retab = [5772 5800 6000 6500 7000  8000 10000 15000 20000 35000 50000 75000 100000]; 
  % first point : from known result
  kinftab = 1.02;
  ksuptab = 1.02; 
  % second point : guess 
  Re = Retab(2)
  kinf = fzero(@(alpha)(real(EV_Poiseuille(Re,alpha))),0.97)
  ksup = fzero(@(alpha)(real(EV_Poiseuille(Re,alpha))),1.05)
  kinftab = [kinftab kinf];
  ksuptab = [ksuptab ksup];
  
  for Re = Retab(3:end)
     Re
     kinf = fzero(@(alpha)(real(EV_Poiseuille(Re,alpha))),kinftab(end))
     ksup = fzero(@(alpha)(real(EV_Poiseuille(Re,alpha))),ksuptab(end))
     kinftab = [kinftab kinf];
     ksuptab = [ksuptab ksup];
  end
  
  figure(3);
  plot(Retab,kinftab,'b',Retab,ksuptab,'b');
 xlabel('Re');ylabel('k');
 set(gcf,'paperpositionmode','auto');
 print('-dpng','-r100','PlanePoiseuille_NeutralCurve.png');

Marginal stability curve

Function performing the eigenvalue computation

function [s,U] = EV_Poiseuille(Re,alpha)

global y dy dyy Z I INT n

dx=i*alpha*I; dxx=-alpha^2*I;
Delta=dxx+dyy;

% base flow
U=1-y.^2; Uy=-2*y;
S=-diag(U)*dx+Delta/Re;

% the matrices
A=[ ...
    S, -diag(Uy), -dx; ...
    Z, S, -dy; ...
    dx, dy, Z];
E=blkdiag(I,I,Z);


% locations on the grid
u=1:n; v=n+1:2*n; p=2*n+1:2*n;

% boundary conditions
III=eye(3*n);
loc=[u(1) u(n) v(1) v(n) ];  
C=III(loc,:);
E(loc,:)=0;  A(loc,:)=C;



% 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)=[];

if(nargout==1) 
    s=s(1); 
    % this is a trick to allow to use the function as s=EV_Poiseuille(Re,alpha). 
    % This allows to pass the function as a handle for fzero.
end

end % function EV_Poiseuille