THE ROTATING POLYGON INSTABILITY OF POTENTIAL VORTEX

This program computes the stability properties of a potential vortex in

The script reproduces results from figures 4, 5 and 6 from Mougel et al. (2018)

Contents

close all;
addpath([fileparts(fileparts(pwd)), '/SOURCES_MATLAB']);
SF_Start('verbosity',2,'ffdatadir','./WORK/');
SF_core_arborescence('cleanall');
SF_core_setopt('storagemode',1);
system('mkdir FIGURES');
figureformat = 'png';
verbosity = 0;

SF_core_setopt('eigensolver','ARPACK'); % this set of programs uses arpack ; to be rewritten someday with Slepc

NOTICE  - Initialization already done

CHAPTER 0 : creation of initial mesh

a=.3;
xi = .25;
rhog = 1;
R = 1;
gamma = 0;
GammaBAR = xi*sqrt(2*rhog*a*R^2/(R^2-xi^2-2*xi^2*log(R/xi)))

g = 1;
H = a;
xi*sqrt(2*g*H*R^2/(R^2-xi^2-2*xi^2*log(R/xi)))


density=60;

freesurf = SF_Init('Mesh_PotentialVortex.edp','Params',[a xi density],'problemtype','3DFreeSurfaceStatic');

GammaBAR =

    0.2215


ans =

    0.2215

CHAPTER 1 : Eigenvalue computation for m=3

warning : in Mougel we have modal dependance with the form exp( i (m theta - omega t) ) and omega is positive (real for pure waves). here the form is exp ( lambda t + i m theta) . So omegar = - lambdai and omegai = lambda_r . This is why the guess is -3i and the eigenvalues are transposed for comparison.

[evm3,emm3] =  SF_Stability(freesurf,'nev',15,'m',3,'shift',-3i);

Plots of the eigenmodes for these parameters

2d plots

figure(10);
subplot(2,2,1);
n=1;
SF_Plot(emm3(n),'phi.im','xlim',[0 1],'ylim',[0 .5]);hold on;SF_Plot_ETA(emm3(n),'Amp',0.05);
title(['Mode G0 : omega = ',num2str(-imag(evm3(n)))])
hold off;
subplot(2,2,2);
n=3;
SF_Plot(emm3(n),'phi.im','xlim',[0 1],'ylim',[0 .5]);hold on;SF_Plot_ETA(emm3(n),'Amp',0.05);
title(['Mode G1 : omega = ',num2str(-imag(evm3(n)))]);xlim([.2 1]);ylim([0 .5]);
hold off;
subplot(2,2,3);
n=2;
SF_Plot(emm3(n),'phi.im','xlim',[0.2 1],'ylim',[0 .5]);hold on;SF_Plot_ETA(emm3(n),'Amp',0.05);
hold on;
title(['Mode C0 : omega = ',num2str(-imag(evm3(n)))]);xlim([.2 1]);ylim([0 .5]);
hold off;
subplot(2,2,4)
n=4;
SF_Plot(emm3(n),'phi.im','xlim',[0 1],'ylim',[0 .5]);hold on;SF_Plot_ETA(emm3(n),'Amp',0.05);
title(['Mode C1 : omega = ',num2str(-imag(evm3(n)))]);xlim([.2 1]);ylim([0 .5]);
hold off;
hold off;
box on; pos = get(gcf,'Position'); pos(3)=pos(4)*1.5;set(gcf,'Position',pos); % resize aspect ratio
saveas(gca,'FIGURES/POLYGONS_modes',figureformat);
pause(0.1);

evm3FIG = 1i*evm3(1:4).' % transpose, not conjugate
figure(22);hold on;
plot(0.25,real(evm3FIG(1)),'ro',0.25,real(evm3FIG(2)),'ro',0.25,real(evm3FIG(3)),'ro',0.25,real(evm3FIG(4)),'ro');


% 3D plots

figure(11);
subplot(2,2,1);SF_Plot_ETA(emm3(1),'Amp',0.05,'dim','3D');title(['Mode G0 : omega = ',num2str(-imag(evm3(1)))]);
subplot(2,2,2);SF_Plot_ETA(emm3(3),'Amp',0.05,'dim','3D');title(['Mode G1 : omega = ',num2str(-imag(evm3(3)))]);
subplot(2,2,3);SF_Plot_ETA(emm3(2),'Amp',0.05,'dim','3D');title(['Mode C0 : omega = ',num2str(-imag(evm3(2)))]);
subplot(2,2,4);SF_Plot_ETA(emm3(4),'Amp',0.05,'dim','3D');title(['Mode C1 : omega = ',num2str(-imag(evm3(4)))]);

pause(0.1);

evm3FIG =

   2.8333 + 0.0000i   3.7210 + 0.0000i   4.2695 + 0.0000i   1.3374 - 0.0000i

CHAPTER 2a : loop for xi = [.25 , .35] by increasing values

freesurf = SF_Init('Mesh_PotentialVortex.edp','Params',[a .25 density],'problemtype','3DFreeSurfaceStatic');
GammaBAR = xi*sqrt(2*rhog*a*R^2/(R^2-xi^2-2*xi^2*log(R/xi)));
evm3 =  SF_Stability(freesurf,'nev',20,'m',3,'shift',-3i); % without cont for initiating branches

tabxi = .25:.0025:.35;
tabEVm3 = [];
figure(20);
title('A few free surface shapes for potential vortex');hold on;

for xi = tabxi
    freesurf = SF_Init('Mesh_PotentialVortex.edp','Params',[a xi density],'problemtype','3DFreeSurfaceStatic');
    GammaBAR = xi*sqrt(2*rhog*a*R^2/(R^2-xi^2-2*xi^2*log(R/xi)));
    figure(20);plot(freesurf.mesh.rsurf,freesurf.mesh.zsurf); hold on;pause(0.1);
    evm3 =  SF_Stability(freesurf,'nev',20,'m',3,'sort','cont','shift',-3i);
    tabEVm3 = [tabEVm3 evm3];
end

% PLOTS
figure(22);hold on;
for num=1:8
    plot(tabxi,-imag(tabEVm3(num,:)),'b-','LineWidth',2);
end
figure(23);hold on;
for num=1:8
    plot(tabxi,real(tabEVm3(num,:)),'b-','LineWidth',2);
end

pause(0.1);

CHAPTER 2b : loop for xi = [.35 , .7] by increasing values

tabxi = .35:.025:.7;
tabEVm3 = [];
for xi = tabxi
    freesurf = SF_Init('Mesh_PotentialVortex.edp','Params',[a xi density],'problemtype','3DFreeSurfaceStatic');
    GammaBAR = xi*sqrt(2*rhog*a*R^2/(R^2-xi^2-2*xi^2*log(R/xi)));
    figure(20);plot(freesurf.mesh.rsurf,freesurf.mesh.zsurf); hold on;pause(0.1);
    evm3 =  SF_Stability(freesurf,'nev',20,'m',3,'sort','cont','shift',-3i);
    tabEVm3 = [tabEVm3 evm3];
end

% PLOTS
figure(22);hold on;
for num=1:8
    plot(tabxi,-imag(tabEVm3(num,:)),'b-','LineWidth',2);
end
figure(23);hold on;
for num=1:8
    plot(tabxi,real(tabEVm3(num,:)),'b-','LineWidth',2);
end

CHAPTER 2c : loop for xi = [.25 , .1] by decreasing values

freesurf = SF_Init('Mesh_PotentialVortex.edp','Params',[a .25 density],'problemtype','3DFreeSurfaceStatic');
evm3 =  SF_Stability(freesurf,'nev',20,'m',3,'shift',-3i); % without cont for initiating branches

tabxi = .25:-.001:.1;
tabEVm3 = [];
for xi = tabxi
    freesurf = SF_Init('Mesh_PotentialVortex.edp','Params',[a xi density],'problemtype','3DFreeSurfaceStatic');
    GammaBAR = xi*sqrt(2*rhog*a*R^2/(R^2-xi^2-2*xi^2*log(R/xi)));
    figure(20);
    plot(freesurf.mesh.rsurf,freesurf.mesh.zsurf); hold on;
    pause(0.1);
    evm3 =  SF_Stability(freesurf,'nev',20,'m',3,'sort','cont','shift',-3i);
    tabEVm3 = [tabEVm3 evm3];
end

% PLOTS
figure(22);hold on;
for num=1:8
    plot(tabxi,-imag(tabEVm3(num,:)),'b','LineWidth',2);
end
title('\omega_r(\xi) for H/R = .3 and m=3');
xlabel('\xi');ylabel('\omega_r');
ylim([0 5]);
box on; pos = get(gcf,'Position'); pos(4)=pos(3)*1;set(gcf,'Position',pos); % resize aspect ratio
set(gca,'FontSize', 14);
saveas(gca,'FIGURES/POLYGONS_omega',figureformat);

figure(23);hold on;
for num=1:8
    plot(tabxi,real(tabEVm3(num,:)),'b-','LineWidth',2);
end
title('\omega_i(\xi) for H/R = .3 and m=3');
xlabel('\xi');ylabel('\omega_r');
ylim([0  .1]);
box on; pos = get(gcf,'Position'); pos(4)=pos(3)*1;set(gcf,'Position',pos); % resize aspect ratio
set(gca,'FontSize', 14);
saveas(gca,'FIGURES/POLYGONS_sigma',figureformat);

% [[PUBLISH]] (this tag is to enable automatic publication as html ; don't touch this line)

%
fclose(fopen('SCRIPT_POLYGONS.success','w+'));

Published with MATLAB® R2020a