sandbox/easystab/kelvin_helmholtz_gerris.m

    Kelvin-Helmholtz : validation with Gerris

    This code computes the instability of the Kelvin-Helmholtz instability of a shear layer, just like kelvin_helmholtz_hermite.m, but here the idea is to generate the velocity field for a nonlinear simulation with Gerris.

    The parameter file for Gerris is kelvin_helmholtz.gfs. We use a domain built out of two square boxes each of size Lx.

    Dependency:

    Files created:

    • eevo showing the evolution in time of the perturbation energy in Gerris
    • u.cgd and v.cgd the initial condition for Gerris
    clear all; clf;
Lx=10;
alpha=2*pi/Lx;    % the wave number
Re=1000;    % the Reynolds number
N=100;      % the number of grid points
L=2*Lx;     % height of domain in gerris

    Here we do a scaling of the y axis and associated differentiation matrices such that the domain in this file corresponds to the height of the domain in Gerris.

    % Hermite differentiation matrices
scale=2*13.4064873381/L;
[y, DM] = herdif(N,2,1);
D=DM(:,:,1)*scale;
DD=DM(:,:,2)*scale^2;
y=y/scale;
Z=zeros(N,N); I=eye(N); 

% renaming the differentiation matrices
dy=D; dyy=DD;
dx=i*alpha*I; dxx=-alpha^2*I;
Delta=dxx+dyy;

% base flow
U=tanh(y); 
Uy=(1-tanh(y).^2);

% the matrices
S=-diag(U)*dx+Delta/Re;
A=[ ...
    S, -diag(Uy), -dx; ...
    Z, S, -dy; ...
    dx, dy, Z];
E=blkdiag(I,I,Z);

% computing eigenmodes 
disp('computing eigenmodes')
[UU,S]=eig(A,E);
s=diag(S);  [t,o]=sort(-real(s)); s=s(o); UU=UU(:,o);
rem=abs(s)>1000; s(rem)=[]; UU(:,rem)=[];

    Initial condition for Gerris

    From now on, we prepare the 2D field for the initial condition of Gerris. Nx is the nimber o grid points in x. And here we select the first eigenmode, which is the least stable one.

    % showing the velocity field
Nx=100;
u=1:N; v=u+N; p=v+N; 
q=UU(:,1); 
Lx=2*pi/alpha;  x=linspace(-Lx/2,Lx/2,Nx);

    We need to translate vertically the axis for Gerris, because we use two boxes and the y=0 is in the center of box one which is the bottom one. This way, the shear layer is located just at the connection between the two boxes of Gerris.

    % translate y axis for gerris
y=y+Lx/2;

    Here we chose the relative amplitude of the eigenmode and the base flow (the base flow has max u velocity equal to 1). Chosing an amplitude too large means that the evolution will be nonlinear right from the beginning, which we don’t want because we would like to compare the evolution of the energy in Gerris and see that as long as the perturbatino is small, it is the same as in this code. If we start on the other hand with an amplitude too small, grid resolution in Gerris and round-off error may disturb the evolution of our eigenmode.

    To expand to physical space, we do like in kelvin_helmholtz_hermite.m#velocity-field.

    % scale mode amplitude 
q=0.01*q/max(abs(q(u))); 

% expand to physical space
qphys=2*real(q*exp(i*alpha*x));

% add the base flow to the perturbations
uu=qphys(u,:)+U*ones(1,Nx);
vv=qphys(v,:);
pp=qphys(p,:);

% show the velocity field
selx=1:5:Nx; sely=1:3:N;
quiver(x(selx),y(sely),uu(sely,selx),vv(sely,selx),'k'); hold on
surf(x,y,pp-10,'facealpha',0.5); shading interp;
axis equal;axis([x(1),x(end),y(1),y(end)]);
xlabel('x'); ylabel('y'); title('Velocity field of the initial condition');

    Saving to disk

    There are several ways to save fields for gerris, here we chose the simplest one, the cgd format (cartesian grid data)

    % save to cartesian grid data format for gerris (.cgd)
disp('saving initial condition')
% for u
fid=fopen('u.cgd','w');
fprintf(fid,'%s\n','2 x y');
fprintf(fid,'%u %u\n',Nx,N);
fprintf(fid,'%f ',x);fprintf(fid,'\n');
fprintf(fid,'%f ',y);fprintf(fid,'\n');
fprintf(fid,'%f\n',uu);  
fclose(fid);
% for v
fid=fopen('v.cgd','w');
fprintf(fid,'%s\n','2 x y');
fprintf(fid,'%u %u\n',Nx,N);
fprintf(fid,'%f ',x);fprintf(fid,'\n');
fprintf(fid,'%f ',y);fprintf(fid,'\n');
fprintf(fid,'%f\n',vv);  
fclose(fid);

    Starting Gerris

    The parameter file for the Gerris simulation is in kelvin_helmholtz.gfs.

    We can start Gerris directly from here. The command to be executed in the shell is stored in the string command. The character “!” at the start of a command tells that the command should be executed in the shell. This command, stored in a string, is then evaluated using the function eval.

    In Gerris it is possible to define some of the parameter values directly from the command line, like we do here to set the value of the boxsize. To set it automatically to the value of Lx in the code, we use sprintf which will insert in the string the value of Lx with two digits after the comma.

    If starting Gerris like this does not work on your system, you can type this in the shell (and change the value 10 if you changed Lx in this code):

    gerris2D -m -Dboxsize=10. kelvin_helmholtz.gfs | gfsview2D v.gfv

    The file v.gfv parameterizes the view of the running dimulation by gfsview. If you don’t have gfsview installed (this is more difficult to install than gerris itself), just remove the line “OutputSimulation” in kelvin_helmholtz.gfs and type in the shell

    gerris2D -m -Dboxsize=10. kelvin_helmholtz.gfs

    The Gerris simulation saves on the disc a file eevo which memorizes the time in the first column and the value of the energy of the perturbation as a second column.

    % start gerris simulation
disp('starting gerris')
command=['!gerris2D -m -Dboxsize=' sprintf('%3.2f',Lx) ' kelvin_helmholtz.gfs | gfsview2D v.gfv'];
disp(command);
eval(command)
    A snapshot from the vorticity of the Gerris simulation

    A snapshot from the vorticity of the Gerris simulation

    Validation

    To validate these computations, we compare the evolution of the energy here and in Gerris. The eigenvalue of the mode we use here is stored in s(1). Since the energy is the square of the mode amplitude, and the mode amplitudes grows in time like \displaystyle \exp(s_rt) where s_r is the real part of the eigenvalue s, then the energy grows like \displaystyle \exp(2s_rt) We scale this energy so that the linear one and the nonlinear one are initialy equal.

    % validation
clf;
a=load('eevo');
t=a(:,1); e=a(:,2);
etheo=exp(2*real(s(1))*t); 
etheo=etheo/etheo(1)*e(1); 

semilogy(t,e,'b.-',t,etheo,'r');
xlabel('time'); ylabel('energy of perturbation')
legend('gerris','easystab');
title('Comparison Gerris/easystab')
grid on

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','kelvin_helmholtz.png')
    Comparison of easystab and Gerris

    Comparison of easystab and Gerris

    Exercices/Contributions

    • Please change this code such that we can have several wavelengthes of the instability in the initial condition, to have a look at the phenomenon of vortex-pairing
    • Please validate the code for several values of the viscosity and of the wavelength
    • Please do these validation with Gerris for other instabilities, like poiseuille_uvp.m or couette_uvwp.m, or pipe_sym.m or even rayleigh_benard.m for which you will need to initiate as well the temperature field.

    %}