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
- Chapter 1 : builds or recovers base flow
- Chapter 2 : plot base flow (something resembling fig. 4)
- CHAPTER 3 : Forced structures for Re = 2000
- 3bb : Plot pressure of the forced structures (re and im)
- 3c : Here is how to plot the VORTICITY (real/imaginary parts, on a more limited range)
- Chapter 4 : compute Eigenmodes
- Chapter 4b : FIGURES FOR EIGENMODES (figure 16 of the paper, p and ux)
- Chapter 4c : FIGURES FOR EIGENMODES (re, Im of p )
- 4d Here is how to plot the VORTICITY (real/imaginary parts, on a more limited range)
- Chapter 5 : Compute Adjoint Eigenmodes / sensitivity (readapt mesh)
- Chapter 5b : FIGURES (fig 17 of the paper)
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 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+'));
