sandbox/easystab/KH_temporal_inviscid.m

    This document belongs to the easystab project, check the main page of the project 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 6

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

    Kelvin-Helmholtz instability of a shear layer

    (temporal analysis, viscous case, discretization based on Rayleigh Equation)

    We study the stability of a parallel shear flow U(y) in the inviscid case. We solve the Rayleigh equation, written as follows :

    \displaystyle (\bar{U} - c) (\partial_y^2 - k^2) \hat \psi - \partial_y^2 \bar{U} \hat \psi = 0

    function [] = main()
close all;
global y dy dyy w Z I U Uy Uyy
loopk = 1; % set to 0 to skipp the loop over k

    Derivation matrices

    Here we wish to solve a problem in an infinite domain. We use Chebyshev discretization with stretching. See differential_equation_infinitedomain.m to see how this works.

    % Numerical parameters
N=201;      % the number of grid points
discretization = 'chebInfAlg'; 
[dy,dyy,w,y] = dif1D(discretization,0,3,N,0.9999);
Z=zeros(N,N); I=eye(N);

    Base flow

    The base flow is U(y) = tanh(y).

    We also need to compute its first and second derivatives:

    U=tanh(y); 
Uy=(1-tanh(y).^2);
Uyy = -2*tanh(y).*(1-tanh(y).^2);

    Eigenvalue computation

    We compute the eigenvalues/eigenmodes with the function KH_inviscid, defined at the end of this program.

    % Physical parameters
alpha=.45;    % the wave number


[c,UU] = KH_inviscid(alpha,N);
omega = c*alpha;

    Plot the spectrum :

    figure(1);
for ind=1:length(omega)
  %%%% plot one eigenvalue
  h=plot(real(omega(ind)),imag(omega(ind)),'*'); hold on
  %%%%  assign the corresponding event on click 
  set(h,'buttondownfcn',{@plotmode,UU(:,ind),omega(ind),alpha});
end
xlabel('\omega_r');
ylabel('\omega_i');
title({['Temporal spectrum for k = ',num2str(alpha)], 'Click on eigenvalues to see the eigenmodes'});
set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','KH_temporal_inviscid_spectrum.png');
    ** Figure : Temporal spectrum of a tanh shear layer

    ** Figure : Temporal spectrum of a tanh shear layer

    Plot the unstable eigenmode :

    plotmode([],[],UU(:,1),omega(1),alpha);

pause(0.1);
set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','KH_temporal_inviscid_mode.png');
    ** Figure : Unstable eigenmmode of the tanh shear layer for k=0.5

    ** Figure : Unstable eigenmmode of the tanh shear layer for k=0.5

    Loop over k to draw the amplification rate as function of the wavenumber

    if loopk
    alphatab = 0:.01:1.;
    lambdatab = [];
    for alpha=alphatab
        [s,~] = KH_inviscid(alpha,N);
        lambdatab = [lambdatab alpha*s(1)];
        figure(4);
        plot(alphatab(1:length(lambdatab)),imag(lambdatab),'b');hold on;
        pause(0.01);
    end
    ylabel('\omega_i');
    xlabel('k');
    title('Growth rate as function of k');
    set(gcf,'paperpositionmode','auto');
    print('-dpng','-r80','KH_temporal_inviscid.png');

end%if loopk

end%function main
    ** Figure : results for the Kelvin-Helmholtz instability (temporal) of a tanh shear layer

    ** Figure : results for the Kelvin-Helmholtz instability (temporal) of a tanh shear layer

    Function KH_inviscid

    Here is the function to compute the eigenvalues/eigenmodes.

    function [s,UU] = KH_inviscid(alpha,N)

global y dy dyy I U Uyy

    The Rayleigh equation is written under the form

    \displaystyle c \left[ \partial_y^2 - k^2 \right] \hat \psi = \left[\bar{U} (\partial_y^2 - k^2) - \partial_y^2 \bar{U} \right] \hat \psi

    which is a generalized eigenvalue problem. After discretizing the operators, is can be written under the matricial form:

    \displaystyle c B X = A X

    where X is the discretized version of \hat \psi.

    for more details on the theory please see Lecture notes for chapter 6

    %differential operators
dx=1i*alpha*I; dxx=-alpha^2*I;
Delta=dxx+dyy;

% the matrices
B=Delta;
A=diag(U)*Delta-diag(Uyy);

%Boundary conditions : we use Dirichlet at the boundaries of the stretched domain
indBC=[1,N];
A(indBC,:)=I(indBC,:);
B(indBC,:)=0;  

% compute the eigenmodes 
[UU,S]=eig(A,B);

% sort the eigenvalues by increasing imaginary part
s=diag(S);  [~,o]=sort(-imag(s)); s=s(o); UU=UU(:,o);
rem=abs(s)>1000; s(rem)=[]; UU(:,rem)=[];

end%function KH_inviscid

    Function plotmode

    This function is used for plotting an eigemode when clicking on the corresponding eigenvalue in the “spectrum” plot.

    function [] = plotmode(~,~,psiM,omega,alpha)
    global y dy dyy
    Yrange = 5;
    figure(2);hold off;
    % plot psi(y)
    u=dy*psiM;
    v=-1i*alpha*psiM;
    vorticity = -alpha^2*psiM-dyy*psiM;
    subplot(1,3,1);hold off;
    plot(real(psiM),y,'r-',imag(psiM),y,'r--');hold on;
    plot(real(u),y,'b-',imag(u),y,'b--');hold on;
    plot(real(v),y,'g-',imag(v),y,'g--');hold on;
    ylabel('y'); ylim([-Yrange,Yrange]);
%legend({'$Re(\hat \psi)$','$Im(\hat \psi)$','$Re(\hat u)$','$Im(\hat u)$','$Re(\hat v)$','$Im(\hat v)$'},'Interpreter','latex');
    legend({'Re(\psi)','Im(\psi)','Re(u)','Im(u)','Re(v)','Im(v)'});
    title('Structure of the eigenmode');
    % plot 2D reconstruction
    Lx=2*pi/alpha; Nx =30;  x=linspace(-Lx/2,Lx/2,Nx);
    psiM(abs(y)>Yrange,:)=[];
    u(abs(y)>Yrange,:)=[];
    v(abs(y)>Yrange,:)=[];
    vorticity(abs(y)>Yrange,:)=[];
    yy = y;
    yy(abs(y)>Yrange)=[];
    psipsi = 2*real(psiM*exp(1i*alpha*x));
    uu=2*real(u*exp(1i*alpha*x));
    vv=2*real(v*exp(1i*alpha*x));
    vorticityvorticity=2*real(vorticity*exp(1i*alpha*x));
    N = length(psiM);
    sely=1:2:N;
    subplot(1,3,2:3); hold off;
    contourf(x,yy,uu,10); hold on; 
    quiver(x,yy(sely),uu(sely,:),vv(sely,:),'k'); hold on;
    %axis equal;
    xlabel('x'); ylabel('y'); title({'Structure of the eigenmode (2D reconstruction)',['for k = ',num2str(alpha) , ' ; omega = ',num2str(omega)]});
end% function plotmode

    Exercices/Contributions

    • Please check the convergence of the results by comparing with other discretization methods (e.g. Hermite)
    • Please modify the program to treat the case of a 2D jet defined as U(y) = 1/cosh(y)^N where N is an integer defining the ‘steepness’ of the profile (try for instance N=1 for a smooth profile and N=10 for a very steep profile). Compare the latter case with theoretical results for a “top-hat” jet. You may need to adapt the mesh parameters to have enough points in the region of the steep gradients !

    %}