Instability and impedance of a jet flowing across a circular aperture in a thick plate (beta = 1)

THIS SCRIPT GENERATES PLOTS FOR THE FORCED STRUCTURES AND THE EIGENMODES FOR THE FLOW THROUGH A HOLE IN A THICK PLATE WITH ASPECT RATIO beta=1

REFEFENCE : Fabre, Longobardi, Citro & Luchini, JFM, 2019

This script reproduces figures 11, 12, 16 and 17 of the paper (and a few more)

You may also have a look at SCRIPT_chi03.m which reproduces figures for the case chi=0.3.

Contents

Chapter 0 : Initialization

addpath([fileparts(fileparts(pwd)), '/SOURCES_MATLAB']);
SF_Start('verbosity',2);
SF_DataBase('open','./WORK_chi1/');
subplot = @(m,n,p) subtightplot (m, n, p, [0.05 0.05], [0.05 0.03], [0.05 0.01]);
mkdir('FIGURES')
if ~SF_core_detectlib('SLEPc')
    error('Error : This case requires SLEPC ! ')
end
NOTICE  - Initialization already done

Chapter 1 : builds or recovers base flow

chi = 1;

bf = SmartMesh_Hole(chi); % or your own sequence for mesh/bf generation
WARNING -  Your program uses depreciated syntax 'Params' to pass parameters 
WARNING -  It is advised to switch to new method using 'Options', and modify your Freefem solvers to detect parameters using 'getARGV'  
WARNING -  Folder ./WORK_chi1/MESHBF not found : creating it !
WARNING -  baseflow.iter = 0 ! it seems your baseflow was projected after mesh adaptation but not recomputed !
WARNING - Using Complex mapping with gamma = 0.3
WARNING - Using Complex mapping with gamma = 0.3
WARNING -  baseflow.iter = 0 ! it seems your baseflow was projected after mesh adaptation but not recomputed !
WARNING - Using Complex mapping with gamma = 0.3
WARNING - Using Complex mapping with gamma = 0.3

Chapter 2 : plot base flow (something resembling fig. 4)

Re = 1600;
bf = SF_BaseFlow(bf,'Re',Re);
SF_Plot(bf,'ux','title','Base Flow','colormap','redblue',...
           'xlim',[-3 3],'ylim',[-1.5 1.5],'contour','on','clevels',[0 0],...
           'boundary','on','bdlabels',2,'bdcolors','k','symmetry','XS');
hold on;
pause(0.1);

% Note that the actual figure 4 of the paper is produced using tecplot.
% To export the data to tecplot format, do the following :
exportFF_tecplot(bf,'BaseFlow.plt');
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
PROGRAM exportFF_tecplot : exporting data in tecplot format
P1 field : vort
P1 field : psi
/bin/bash: tec360 : commande introuvable
END PROGRAM exportFF_tecplot

CHAPTER 3 : Forced structures for Re = 2000

Computes forced structures for 5 values of omega
 C1= 0,95   S1=1,85      H2=2,5      C2=2,7      S2=3,8
bf = SF_BaseFlow(bf,'Re',2000);
Params = [5 17 1.25 .5 5 17]; % choice of parameters to use complex mapping + stretching
% Params = [5 1e30 1.25 .5 20 1e30];% alternative choice to use complex mapping without stretching

omega = 0.95;
foA = SF_LinearForced(bf,omega,'mappingdef','jet','mappingparams',Params);

omega = 1.85;
foB = SF_LinearForced(bf,omega,'mappingdef','jet','mappingparams',Params);

omega = 2.5;
foC = SF_LinearForced(bf,omega,'mappingdef','jet','mappingparams',Params);

omega = 2.7;
foD = SF_LinearForced(bf,omega,'mappingdef','jet','mappingparams',Params);

omega = 3.8;
foE = SF_LinearForced(bf,omega,'mappingdef','jet','mappingparams',Params);
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5

3bb : Plot pressure of the forced structures (re and im)

(figure 11 of the paper)
figure;
Re = 2000;

subplot(5,2,1);
SF_Plot(foA,'p','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',1);
           text(-2.8,3.2,'(C1)');
subplot(5,2,2);
SF_Plot(foA,'p.im','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',1);
subplot(5,2,3);
SF_Plot(foB,'p','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',.3);
           text(-2.8,3.2,'(S1)');
subplot(5,2,4);
SF_Plot(foB,'p.im','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',.3);
subplot(5,2,5);
SF_Plot(foC,'p','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',.3);
           text(-2.8,3.2,'(H2)');
subplot(5,2,6);
SF_Plot(foC,'p.im','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',.3);
subplot(5,2,7);
SF_Plot(foD,'p','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',.3);
           text(-2.8,3.2,'(C2)');
subplot(5,2,8);
SF_Plot(foD,'p.im','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',.3);
subplot(5,2,9);
SF_Plot(foE,'p','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',.3);
           text(-2.8,3.2,'(S2)');
subplot(5,2,10);
SF_Plot(foE,'p.im','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',.3);

pos = get(gcf,'Position'); pos(3) = 800; pos(4)=pos(3)*1.18;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/ForcedModes_P_chi', num2str(chi), '_Re',num2str(Re),'.png'],'png');
saveas(gcf,['FIGURES/ForcedModes_P_chi', num2str(chi), '_Re',num2str(Re), '.fig'],'fig');
saveas(gcf,['FIGURES/ForcedModes_P_chi', num2str(chi), '_Re',num2str(Re), '.tif'],'tif');

3c : Here is how to plot the VORTICITY (real/imaginary parts, on a more limited range)

(Figure 12 of the paper)

figure;
subplot = @(m,n,p) subtightplot (m, n, p, [0.05 0.05], [0.05 0.03], [0.05 0.01]);

subplot(5,2,1);
SF_Plot(foA,'vort.re','boundary','on','colormap','redblue','colorrange',[-50 50],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.1,1.45,'(C1)');
subplot(5,2,2);
SF_Plot(foA,'vort.im','boundary','on','colormap','redblue','colorrange',[-50 50],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,i}');
subplot(5,2,3);
SF_Plot(foB,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.1,1.45,'(S1)');
subplot(5,2,4);
SF_Plot(foB,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
subplot(5,2,5);
SF_Plot(foC,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.1,1.45,'(H2)');
subplot(5,2,6);
SF_Plot(foC,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,i}');
subplot(5,2,7);
SF_Plot(foD,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.1,1.45,'(C2)');
subplot(5,2,8);
SF_Plot(foD,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,i}');
subplot(5,2,9);
SF_Plot(foE,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.1,1.45,'(S2)');
subplot(5,2,10);
SF_Plot(foE,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,i}');


pos = get(gcf,'Position'); pos(3) = 800; pos(4)=pos(3)*1.18;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re),'.png'],'png');
saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re), '.fig'],'fig');
saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re), '.tif'],'tif');

Chapter 4 : compute Eigenmodes

Re =1500;
Params = [5 17 1.25 .5 5 17];
bf = SF_BaseFlow(bf,'Re',Re);
[ev,em2D] = SF_Stability(bf,'solver','StabAxi.edp','shift',-2.1i,'m',0,'nev',1,'type','D',...
               'mappingdef','jet','mappingparams',Params);

[ev,em2A] = SF_Stability(bf,'shift',-2.1i,'m',0,'nev',1,'type','A',...
               'mappingdef','jet','mappingparams',Params);

em2S = SF_Sensitivity(bf,em2D,em2A);

Re =1700;
bf = SF_BaseFlow(bf,'Re',Re);
[ev,em3D] = SF_Stability(bf,'shift',-4.1i,'m',0,'nev',1,'type','D',...
               'mappingdef','jet','mappingparams',Params);

[ev,em3A] = SF_Stability(bf,'shift',-4.14i,'m',0,'nev',1,'type','A',...
               'mappingdef','jet','mappingparams',Params);

em3S = SF_Sensitivity(bf,em3D,em3A);
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING -  Sensitivity : here this is not normalized ; amplitude is arbitrary
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING -  Sensitivity : here this is not normalized ; amplitude is arbitrary

Chapter 4b : FIGURES FOR EIGENMODES (figure 16 of the paper, p and ux)

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.05 0.05], [0.05 0.01], [0.05 0.01]);
if ~make_it_tight,  clear subplot;  end

figure(33);
subplot(2,2,1);
SF_Plot(em2D(1),'p','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',1);
text(-2.8,3.2,'(a)');
subplot(2,2,2);hold on;
SF_Plot(em2D(1),'vort','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_{x,r}','logsat',1);
text(-2.8,3.2,'(b)');
subplot(2,2,3);
SF_Plot(em3D,'p','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',1);
box on; text(-2.8,3.2,'(c)');
subplot(2,2,4);
SF_Plot(em3D,'vort','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_{x,r}','logsat',1);
box on; text(-2.8,3.2,'(d)');

pos = get(gcf,'Position'); pos(3) = 800;pos(4)=pos(3)*.4;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/EigenModes_PU_chi', num2str(chi),'.png'],'png');
saveas(gcf,['FIGURES/EigenModes_PU_chi', num2str(chi),'.fig'],'fig');
saveas(gcf,['FIGURES/EigenModes_PU_chi', num2str(chi),'.tif'],'tif')
pause(0.1);

Chapter 4c : FIGURES FOR EIGENMODES (re, Im of p )

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.05 0.05], [0.05 0.01], [0.05 0.01]);
if ~make_it_tight,  clear subplot;  end

figure(35);
subplot(2,2,1);
SF_Plot(em2D(1),'p','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',1);
text(-2.8,3.2,'(a)');
subplot(2,2,2);hold on;
SF_Plot(em2D(1),'p.im','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',1);
text(-2.8,3.2,'(b)');
subplot(2,2,3);
SF_Plot(em3D,'p','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',1);
box on; text(-2.8,3.2,'(c)');
subplot(2,2,4);
SF_Plot(em3D,'p.im','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',1);
box on; text(-2.8,3.2,'(d)');

pos = get(gcf,'Position'); pos(3) = 800;pos(4)=pos(3)*.4;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/EigenModes_P_chi', num2str(chi),'.png'],'png');
saveas(gcf,['FIGURES/EigenModes_P_chi', num2str(chi),'.fig'],'fig');
saveas(gcf,['FIGURES/EigenModes_P_chi', num2str(chi),'.tif'],'tif');
pause(0.1);

4d Here is how to plot the VORTICITY (real/imaginary parts, on a more limited range)

figure;
subplot = @(m,n,p) subtightplot (m, n, p, [0.1 0.1], [0.05 0.05], [0.05 0.05]);

subplot(2,2,1);
SF_Plot(em2D,'vort.re','boundary','on','colormap','redblue','colorrange',[-50 50],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.2,1.5,'(a)');
subplot(2,2,2);
SF_Plot(em2D,'vort.im','boundary','on','colormap','redblue','colorrange',[-50 50],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,i}');
             text(-2.2,1.5,'(b)');
subplot(2,2,3);
SF_Plot(em3D,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.8,3.5,'(b)');
subplot(2,2,4);
SF_Plot(em3D,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-2.2 0.2],'ylim',[.6 1.4],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{z,r}');
             text(-2.8,3.5,'(b)');

pos = get(gcf,'Position'); pos(3) = 800;pos(4)=pos(3)*.4;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re),'.png'],'png');
saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re), '.fig'],'fig');
saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re), '.tif'],'tif');

if 0

Chapter 5 : Compute Adjoint Eigenmodes / sensitivity (readapt mesh)

Re =1500;
bf = SF_BaseFlow(bf,'Re',Re);
%[ev,em2S,em2D,em2A] = SF_Stability(bf,'shift',-2.1i,'m',0,'nev',1,'type','S',...
%               'mappingdef','jet','mappingparams',Params);
[ev,em2D] = SF_Stability(bf,'shift',-2.1i,'m',0,'nev',1,'type','D',...
               'mappingdef','jet','mappingparams',Params);

[ev,em2A] = SF_Stability(bf,'shift',-2.1i,'m',0,'nev',1,'type','A',...
               'mappingdef','jet','mappingparams',Params);

em2S = SF_Sensitivity(bf,em2D,em2A);

Re =1700;
bf = SF_BaseFlow(bf,'Re',Re);

[ev,em3D] = SF_Stability(bf,'shift',-4.14i,'m',0,'nev',1,'type','D',...
               'mappingdef','jet','mappingparams',Params);

[ev,em3A] = SF_Stability(bf,'shift',-4.14i,'m',0,'nev',1,'type','A',...
               'mappingdef','jet','mappingparams',Params);

em3S = SF_Sensitivity(bf,em3D,em3A);
end

Chapter 5b : FIGURES (fig 17 of the paper)

figure(34);
subplot(2,2,1);
SF_Plot(em2A,'ux','colormap','redblue','colorrange',[-.5 .5],'xlim',[-3.5 .5],'ylim',[0.5 2],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_{x,adj,r}');
text(-3.45,2.1,'(a)');
subplot(2,2,2);
SF_Plot(em2S,'sensitivity','colormap','ice','colorrange',[0 0.04],'xlim',[-2.25 1.75],'ylim',[0 1.5],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','S_w','colorbar','eastoutside');
text(-2.2,1.6,'(b)');
subplot(2,2,3);
SF_Plot(em3A,'ux','colormap','redblue','colorrange',[-.2 .2],'xlim',[-3.5 .5],'ylim',[0.5 2],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_{x,adj,r}');
box on; text(-3.45,2.1,'(c)');
subplot(2,2,4);
SF_Plot(em3S,'sensitivity','colormap','ice','colorrange',[0 0.04],'xlim',[-2.25 1.75],'ylim',[0 1.5],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','S_w');
box on; text(-2.2,1.6,'(d)');
pos = get(gcf,'Position'); pos(3) = 800;pos(4)=pos(3)*.4;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/EigenModes_Adj_chi', num2str(chi),'.png'],'png');
saveas(gcf,['FIGURES/EigenModes_Adj_chi', num2str(chi),'.fig'],'fig');
saveas(gcf,['FIGURES/EigenModes_Adj_chi', num2str(chi),'.tif'],'tif');
pause(0.1);


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

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