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+'));


