Harmonically forced flow (and impedance) for jet through a hole with aspect ratio 0.3

This script reproduces part of the results of the paper Fabre, Longobardi, Citro & Luchini (JFM 2020) * Structure of base flow for Re = 1600 (figure 5c from paper) * Forced structures correrponding to frequencies C1, C2 and anti-S1 (figures 8-9 from paper) * Impedance curves for Re = 1600 (figure 7 e-f from paper)

NB this script is in two parts ; the first one performs the computation, the second ones does the post-processing by reading results from the database. if you have previously run and only want to generate the figures, can run only the second part.

Contents

Part I : computations

Initialization

addpath('../../SOURCES_MATLAB');
SF_Start('verbosity',4);
SF_DataBase('open','./WORK_chi03/');
mkdir('FIGURES');
NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead
Warning: Directory already exists. 

Chapter I-1 : builds an adapted mesh and an initial BF

chi = 0.3;
bf = SmartMesh_Hole(0.3); % generate an adapted mesh for Re=1500
Re = 1600;
bf = SF_BaseFlow(bf,'Re',1600,'solver','Newton_Axi.edp');
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_chi03/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 I-2 : Forced structures

Computes forced structures for 3 values of omega

Params = [5 17 1.25 .5 5 17]; % this choice of parameters to plot structures
omega = 2.6;
foA = SF_LinearForced(bf,omega,'solver','LinearForcedAxi_COMPLEX_m0.edp','mappingdef','jet','mappingparams',Params);

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

omega = 8.25;
foC = 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

Chapter I-3 : Impedances curves for this case

Params = [0 17 2.5 .5 5 17]; % this choice of parameters to compute impedance
fo = SF_LinearForced(bf,'omega',[0:.1:8],'mappingdef','jet','mappingparams',Params,'plot','yes');
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5
WARNING -  Folder ./WORK_chi03/IMPEDANCES not found : creating it !

Part II : postprocessing

New initialization in case you execute the script from here:

addpath('../../SOURCES_MATLAB');
SF_Start('verbosity',4);
SF_DataBase('open','./WORK_chi03/');
chi = 0.3;
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Chapter II-1 : Plot BF

bf = SF_Load('BASEFLOWS','last'); % load from database
figure;
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');
pause(0.1);
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

Chapter II-2 : plot forced structures with various representations

first Reload forced structures from database

foA = SF_Load('FORCEDFLOWS',3);
foB = SF_Load('FORCEDFLOWS',4);
foC = SF_Load('FORCEDFLOWS',5);

Plot pressure and axial velocity

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

% a and b
subplot(3,2,1);
SF_Plot(foA,'p','boundary','on','colormap','redblue','colorrange',[-3,3],'xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''','logsat',1);
           text(-2.8,3.5,'(a)');
subplot(3,2,2);
SF_Plot(foA,'ux','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_x','logsat',1);
             text(-2.8,3.5,'(b)');

% c and d
subplot(3,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''','logsat',.5);
           text(-2.8,3.5,'(c)');
subplot(3,2,4);
SF_Plot(foB,'ux','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_x','logsat',.2);
             text(-2.8,3.5,'(d)');

% e and f
subplot(3,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''','logsat',1);
           text(-2.8,3.5,'(e)');
subplot(3,2,6);
SF_Plot(foC,'ux','boundary','on','colormap','redblue','colorrange','cropcentered','xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','u''_x','logsat',1);
             text(-2.8,3.5,'(f)');


 pos = get(gcf,'Position'); pos(4)=pos(3)*.6;set(gcf,'Position',pos);
saveas(gcf,['FIGURES/ForcedModes_PU_chi', num2str(chi), '_Re',num2str(Re),'.png'],'png');
saveas(gcf,['FIGURES/ForcedModes_PU_chi', num2str(chi), '_Re',num2str(Re), '.fig'],'fig');
%saveas(gcf,['FIGURES/ForcedModes_chi', num2str(chi), '_Re',num2str(Re), '.svg'],'svg');

Here is how to plot the PRESSURE

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

% a and b
subplot(3,2,1);
SF_Plot(foA,'p','boundary','on','colormap','redblue','colorrange',[-3,3],'xlim',[-3 5],'ylim',[0 3],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_r','logsat',1);
           text(-2.8,3.5,'(C1)');
subplot(3,2,2);
SF_Plot(foA,'p.im','boundary','on','colormap','redblue','colorrange',[-3,3],'xlim',[-3 5],'ylim',[0 3],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','p''_i','logsat',1);

% c and d
subplot(3,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',.2);
           text(-2.8,3.5,'(S1)');
subplot(3,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',.2);

% e and f
subplot(3,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',1);
           text(-2.8,3.5,'(C2)');
subplot(3,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',1);


 pos = get(gcf,'Position'); pos(4)=pos(3)*.6;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_chi', num2str(chi), '_Re',num2str(Re), '.svg'],'svg');

Plot the pressure and axial velocity along axis

Xa = [-4*chi :.01*chi :2*chi];
figure;

pA = SF_ExtractData(foA,'p',Xa,0);
uA = SF_ExtractData(foA,'ux',Xa,0);
subplot(3,1,1);
plot(Xa,real(uA),'r',Xa,imag(uA),'r--',Xa,real(pA),'b',Xa,imag(pA),'b--');
hold on;
plot([-2*chi -2*chi], [-2 2], 'k:',[0 0], [-2 2], 'k:');
xlabel('x');
legend('u_{x,r}','u_{x,i}','p_r','p_i');

pB = SF_ExtractData(foB,'p',Xa,0);
uB = SF_ExtractData(foB,'ux',Xa,0);
subplot(3,1,2);
plot(Xa,real(uB),'r',Xa,imag(uB),'r--',Xa,real(pB),'b',Xa,imag(pB),'b--');
hold on;
plot([-2*chi -2*chi], [-2 2], 'k:',[0 0], [-2 2], 'k:');
xlabel('x');
legend('u_{x,r}','u_{x,i}','p_r','p_i');

pC = SF_ExtractData(foC,'p',Xa,0);
uC = SF_ExtractData(foC,'ux',Xa,0);
subplot(3,1,3);
plot(Xa,real(uC),'r',Xa,imag(uC),'r--',Xa,real(pC),'b',Xa,imag(pC),'b--');
hold on;
plot([-2*chi -2*chi], [-2 2], 'k:',[0 0], [-2 2], 'k:');
xlabel('x');
legend('u_{x,r}','u_{x,i}','p_r','p_i');

saveas(gcf,['FIGURES/ForcedModes_chi', num2str(chi), '_Re',num2str(Re),'.png'],'png');
saveas(gcf,['FIGURES/ForcedModes_chi', num2str(chi), '_Re',num2str(Re), '.fig'],'fig');
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
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

Here is how to plot the VORTICITY

figure(12);

% a and b
subplot(3,2,1);
SF_Plot(foA,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,r}');
           text(-2.8,3.5,'(a)');
subplot(3,2,2);
SF_Plot(foA,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,i}');
             text(-2.8,3.5,'(b)');

% c and d
subplot(3,2,3);
SF_Plot(foB,'vort','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,r}');
           text(-2.8,3.5,'(c)');
subplot(3,2,4);
SF_Plot(foB,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,i}');
             text(-2.8,3.5,'(d)');

% e and f
subplot(3,2,5);
SF_Plot(foC,'vort','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
               'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,r}');
           text(-2.8,3.5,'(e)');
subplot(3,2,6);
SF_Plot(foC,'vort.im','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,i}');
             text(-2.8,3.5,'(f)');


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

Here is how to plot the VORTICITY (only real part)

figure;
subplot = @(m,n,p) subtightplot (m, n, p, [0.1 0.1], [0.05 0.05], [0.05 0.05]);
% a and b
subplot(3,1,1);
SF_Plot(foA,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,r}');
             text(-.85,1.2,'(C1)');

% c and d
subplot(3,1,2);
SF_Plot(foB,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,r}');
             text(-.85,1.2,'(S1)');

% e and f
subplot(3,1,3);
SF_Plot(foC,'vort.re','boundary','on','colormap','redblue','colorrange',[-100 100],'xlim',[-.9 .3],'ylim',[.675 1.125],...
                'boundary','on','bdlabels',2,'bdcolors','k','cbtitle','\omega''_{\phi,r}');
             text(-.85,1.2,'(C2)');


pos = get(gcf,'Position'); pos(4)=pos(3)*.6*1.2; pos(3) = pos(3)/2*1.2; set(gcf,'Position',pos);
saveas(gcf,['FIGURES/ForcedModes_VORTre_chi', num2str(chi), '_Re',num2str(Re),'.png'],'png');
saveas(gcf,['FIGURES/ForcedModes_VORTre_chi', num2str(chi), '_Re',num2str(Re), '.fig'],'fig')
%saveas(gcf,['FIGURES/ForcedModes_VORT_chi', num2str(chi), '_Re',num2str(Re), '.svg'],'svg')

Chapter II-3 : Plot impedances curves and Nyquist

First reload from database :

fo = SF_LinearForced(bf,'omega',[0:.1:8],'mappingdef','jet','mappingparams',Params);
WARNING - Using Complex mapping with gamma = 0.5
WARNING - Using Complex mapping with gamma = 0.5

plotting impedance curves

fo = SF_Load('IMPEDANCES','last');

clear subplot;
figure; subplot(1,2,1);theplotsymbol='r';
 plot(fo.omega,real(fo.Z),[theplotsymbol,'-'],fo.omega,-imag(fo.Z)./fo.omega,[theplotsymbol,'--']);hold on;
 plot(fo.omega,0*real(fo.Z),'k:','LineWidth',1)
 xlabel('\omega');ylabel('Z_r, -Z_i/\omega');
 title(['Impedance for Re = ',num2str(bf.Re)] );
 subplot(1,2,2);hold on;
 plot(real(fo.Z),imag(fo.Z),[theplotsymbol,'-']); title(['Nyquist diagram for Re = ',num2str(bf.Re)] );
 xlabel('Z_r');ylabel('Z_i');ylim([-10 2]);
 box on; pos = get(gcf,'Position'); pos(4)=pos(3)*.5;set(gcf,'Position',pos);

[[PUBLISH]]

fclose(fopen('SCRIPT_chi03.success','w+'));