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)
- (numerous extra figures !)
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
- Chapter I-1 : builds an adapted mesh and an initial BF
- Chapter I-2 : Forced structures
- Chapter I-3 : Impedances curves for this case
- Part II : postprocessing
- Chapter II-1 : Plot BF
- Chapter II-2 : plot forced structures with various representations
- Plot pressure and axial velocity
- Here is how to plot the PRESSURE
- Plot the pressure and axial velocity along axis
- Here is how to plot the VORTICITY
- Here is how to plot the VORTICITY (only real part)
- Chapter II-3 : Plot impedances curves and Nyquist
- plotting impedance curves
- [[PUBLISH]]
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
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+'));