sandbox/easystab/david/wave_like_Psi_noslip.m
Wave-like pertubations WITH No-Slip CONDITIONS
This code is adapted from wave_like.m from the easystab project.
Modified by D. Fabre to experiment the following things :
Psi formulation compared to u-v-p formulation,
noslip conditions versus no-slip conditions.
Main program :
clear all; clf;
% parameters
N=50; % number of gridpoints
L=1; % Fluid height in y
rho=1; % fluid density
mu=1; % fuid viscosity
alpha=4 ; % wavenumber in x
g=1; % gravity
[D,DD,wy,y]=dif1D('cheb',0,L,N,3);
u=1:N; v=u+N; p=v+N;
% validation
figure(1);
alphavec=linspace(0.01,10,100);
xi = @(lambda,alpha)(sqrt(-alpha^2-lambda));
Fs = @(lambda,alpha)(cosh(alpha*.5)*(-xi(lambda,alpha)*sin(xi(lambda,alpha)*.5)))-alpha*sinh(alpha*.5)*cos(xi(lambda,alpha)*.5);
Fa = @(lambda,alpha)(sinh(alpha*.5)*(xi(lambda,alpha)*cos(xi(lambda,alpha)*.5)))-alpha*cosh(alpha*.5)*sin(xi(lambda,alpha)*.5);
stheos(1)=fzero(@(lambda)(Fs(lambda,alphavec(1))),-40);
stheoa(1)=fzero(@(lambda)(Fa(lambda,alphavec(1))),-80);
for ind=2:length(alphavec)
stheos(ind)=fzero(@(lambda)(Fs(lambda,alphavec(ind))),stheos(ind-1));
stheoa(ind)=fzero(@(lambda)(Fa(lambda,alphavec(ind))),stheoa(ind-1))
end
plot(alphavec,stheos,'r-','DisplayName', 'theory'); hold on
plot(alphavec,stheoa,'r-','DisplayName', 'theory'); hold on
for alpha = 1:10
[s,U] = Wavelike_UVP_noslip(mu,alpha,L,rho,N);
[spsi,Upsi] = Wavelike_psi_noslip(mu,alpha,L,rho,N);
plot(alpha,s(1:4),'b*','DisplayName', 'UVP');
plot(alpha,spsi(1:4),'go','DisplayName', 'Psi')
end
xlabel('alpha');ylabel('exponential growth rate'); title('validation')
hr = findobj('Color','r');
hb = findobj('Color','b');
hg = findobj('Color','g');
hv = [hr(1) hb(1) hg(1)];
legend(hv);
%legend('analytical','UVP','psi');
grid on
set(gcf,'paperpositionmode','auto');
print('-dpng','-r100','wave_like_UVP_Psi_noslip.png');
% show the velocity field of the first three eigenvectors
figure(2);
Lx=2*pi/alpha;
x=linspace(0,Lx,20);
for ind=1:4
subplot(2,2,ind);
q=U(:,ind);
qphys=2*real(q*exp(i*alpha*x));
quiver(x,y(1:2:N),qphys(u(1:2:N),:),qphys(v(1:2:N),:));
axis equal; axis([0,Lx,0,L]);xlabel('x');ylabel('y');
title(['mode' num2str(ind)]);
end
set(gcf,'paperpositionmode','auto');
print('-dpng','-r100','wave_like_Modes_noslip.png');
Here are the eigenvalues, compared to the theory.
On the figure, we show the velocity field for the four least stable eigenmodes. All off course have the same wavelength 2\pi/\alpha in x, but we see that as we go to more decaying modes, they have more and more oscillations in the wall-normal direction. The first mode corresponds to two counter-rotating vortices occupying the whole height between the walls. The second mode is four counter-rotating vortices, since now from y=0 to L, space is made for two vortices, and so on for the following modes.
Function for U-V-P formulation :
function [s,U] = Wavelike_UVP_noslip(mu,alpha,L,rho,N);
% 1D differentiation matrices
[D,DD,wy,y]=dif1D('cheb',0,L,N,3);
I=eye(N); Z=zeros(N,N);
% renaming the matrices
dy=D; dyy=DD;
dx=i*alpha*I; dxx=-alpha^2*I;
Delta=dxx+dyy;
% location vectors
u=1:N; v=u+N; p=v+N;
% System matrices
A=[mu*Delta, Z, -dx; ...
Z, mu*Delta, -dy; ...
dx, dy, Z];
E=[rho*I, Z, Z; ...
Z, rho*I, Z; ...
Z, Z, Z];
% boundary conditions
loc=[u(1),u(N),v(1),v(N)];
II=eye(3*N); DD=blkdiag(dy,dy,dy);
C=[II([u(1),u(N)],:); II([v(1),v(N)],:)];
E(loc,:)=0;
A(loc,:)=C;
% compute 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)=[];
end
Function for Psi formulation (Orr-Sommerfeld) :
function [s,U] = Wavelike_psi_noslip(mu,alpha,L,rho,N);
% differentiation matrices
[y,DM] = chebdif(N,4);
y = L*y/2;
d.y=2/L*DM(:,:,1);
d.yy=(2/L)^2*DM(:,:,2);
d.yyyy=(2/L)^4*DM(:,:,4);
I=eye(N);Z=zeros(N,N);
d.x=i*alpha*I;
% base flow and derivatives
u0=0;
up=0;
upp=0;
% laplacian
k2=alpha^2;
lap=(d.yy-k2*I);
lap2=(d.yyyy-2*k2*d.yy+k2^2*I);
% system matrices
LOS=mu*lap2;
A=[LOS];
E=[lap];
% boundary conditions
psi=1:N;% location vectors
loc=[1,2,N-1,N];
C=[I([1,N],:) ; d.y([1,N],:) ];
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)=[];
end