Acoustic field in a pipe with harmonic forcing at the bottom

This scripts demonstrates the efficiency of StabFem for a linear acoustics problem.

Problem : find the velocity potential such as :

• • along • Sommerfeld radiation condition on ( is the acoustic wavenuber)

Variational formulation : Resolution is made by the FreeFem++ solver LinearForcedAcoustic.edp which is interfaced by the driver SF_LinearForced.m

Additional results on this situation can be found in two more advanced examples available in the StabFem project :

- SCRIPT_PIPE_Comparison_theory.m for comparison of impedances with theory

- SCRIPT_ACOUSTIC_PML_CM.m for demonstration of how to use PML and CM methods for nonreflective boundary conditions

initialisation

close all
SF_Start('verbosity',2);

set(groot, 'defaultAxesTickLabelInterpreter','latex');  % options for figures
set(groot, 'defaultLegendInterpreter','latex'); %idem
set(0,'defaultAxesFontSize',20); % idem
mkdir('html');

Chapter 1 : building an adapted mesh

(See explanations in next chapter about usage of 'Harmonic_Acoustic.edp')

ffmeshInit = SF_Mesh('Mesh_1.edp','Params',10);
Forced = SF_Launch('Harmonic_Acoustic.edp','Options',{'omega',1},'mesh',ffmeshInit,'DataFile','ForcedFlow.txt','Store','FORCEDFLOWS');
Forced1 = SF_Launch('Harmonic_Acoustic.edp','Options',{'omega',.5},'mesh',ffmesh,'DataFile','ForcedFlow.txt','Store','FORCEDFLOWS');
Forced2 = SF_Launch('Harmonic_Acoustic.edp','Options',{'omega',2},'mesh',ffmesh,'DataFile','ForcedFlow.txt','Store','FORCEDFLOWS');
%plot the mesh :

figure;  SF_Plot(ffmeshInit,'mesh','symmetry','ym','boundary','on');
hold on; SF_Plot(ffmesh,'mesh','title','Mesh : Initial (left) and Adapted (right)','boundary','on'); Chapter 2 : Compute and plot the pressure fied with harmonic forcing at the bottom of the tube

omega = 1;
Forced = SF_Launch('Harmonic_Acoustic.edp','Options',{'omega',omega},'mesh',ffmesh,'DataFile','ForcedFlow.txt','Store','FORCEDFLOWS');

% We use the program Harmonic_Acoustic.edp, which computes the solution of the Helmholtz equation
% with forcing at the bottom of the pipe, and writes the results in a file 'ForcedFlow.txt'.
% The program is interfaced using the generic driver SF_Launch, the result is a flowfield object
% containing the forced structure.

% NB : alternatively you can use the driver SF_Linear forced whose syntax is as follows:
% > Forced = SF_LinearForced(ffmesh,'omega',omega,'solver','LinearForcedAcoustic.edp')
% This syntax is depreciated, and it is now recommended to use SF_Launch for single-omega computations.

figure();
SF_Plot(Forced,'p.im','boundary','on','colormap','redblue','cbtitle','Re(p'')',...
'YLim',[-10,20]);
hold on;
SF_Plot(Forced,'p.re','boundary','on','colormap','redblue','symmetry','YM','cbtitle','Im(p'') and Re(p'')',...
'YLim',[-10,20]); Create a movie (animated gif) from this field

h = figure;
filename = 'html/SCRIPT_DEMO_ACOUSTIQUE_movie.gif';
SF_Plot(Forced,'p','boundary','on','colormap','redblue','colorrange',[-1 1],...
'symmetry','YS','cbtitle','p''','colorbar','eastoutside','bdlabels',[1 2 ],'bdcolors','k','Amp',1);
set(gca,'nextplot','replacechildren');
for k = 1:20
Amp = exp(-2*pi*1i*k/20);
SF_Plot(Forced,'p','boundary','on','contour','on','clevels',[-2 :.5 :2],...
'colormap','redblue','colorrange',[-1 1],...
'symmetry','YS','cbtitle','p''','colorbar','eastoutside','bdlabels',[1 2 ],'bdcolors','k','Amp',Amp);
frame = getframe(h);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if k == 1
imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
else
imwrite(imind,cm,filename,'gif','WriteMode','append');
end
end Here is the movie Extract p and u along the symmetry axis

Xaxis = [-10 :.1 :10];
Uyaxis = SF_ExtractData(Forced,'uz',0,Xaxis);
Paxis = SF_ExtractData(Forced,'p',0,Xaxis);
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

Plot p and u along the symmetry axis

figure();
plot(Xaxis,real(Uyaxis),Xaxis,imag(Uyaxis)); hold on;plot(Xaxis,real(Paxis),Xaxis,imag(Paxis));
xlabel('x');
legend('Re(u''z)','Im(u''z)','Re(p'')','Im(p'')');
pause(0.1); Chapter 3 : loop over k to compute the impedance We use the same program "Harmonic_Acoustic.edp" but with a different set of options. In this case the 'DataFile' generated at the end is called "LinearForcedStatistics.ff2m". It is imported and returned as an output a structure 'IMP' containing arrays for omega and Z (and possibly other customizable fields ; see reference manual)

IMP = SF_Launch('Harmonic_Acoustic.edp','Options',{'omegamin',0,'omegamax',2,'Nomega',401},'mesh',ffmesh,'DataFile','LinearForcedStatistics.ff2m','Store','IMPEDANCES')
WARNING -  Folder ./WORK//IMPEDANCES mentioned for database storing does not exist ; creating it

IMP =

struct with fields:

mesh: [1x1 struct]
filename: './WORK/IMPEDANCES/FFDATA_000001.txt'
DataDescription: 'Impedance of a axisymmetric acoustic flow'
datatype: 'ForcedLinear'
datastoragemode: 'columns'
meshfilename: './WORK/MESHES/FFMESH_000005.msh'
ind: [401x1 double]
omega: [401x1 double]
Z: [401x1 double]
R: [401x1 double]

Nb : alternatively you can use the driver "SF_LinearForced" to do the same thing, with a different syntax : IMP = SF_LinearForced(ffmesh,'solver','LinearForcedAcoustic.edp','omega',[0.005:.005:2],'plot','no')

Plot figure;
k=1;
plot(IMP.omega,real(IMP.Z),'b',IMP.omega,imag(IMP.Z),'b--');
hold on;

title(['Impedance $Z_r$ and $Z_i$'],'Interpreter','latex','FontSize', 30)
xlabel('$k$','Interpreter','latex','FontSize', 30);
ylabel('$Z_r,Z_i$','Interpreter','latex','FontSize', 30);
set(findall(gca, 'Type', 'Line'),'LineWidth',2);
pause(0.1); plot in semilog

figure;
plot(IMP.omega,abs(IMP.Z),'k--');
xlabel('b'); ylabel('|Z|');
xlabel('$k$','Interpreter','latex','FontSize', 30);
ylabel('$|Z|$','Interpreter','latex','FontSize', 30);
title(['Impedance $|Z|$'],'Interpreter','latex','FontSize', 30)
set(findall(gca, 'Type', 'Line'),'LineWidth',2);

pause(0.1); plot reflection coefficient

figure;
semilogy(IMP.omega,IMP.R,'b--');
xlabel('$k$','Interpreter','latex','FontSize', 30);
ylabel('$R_i$','Interpreter','latex','FontSize', 30);
title(['Reflection coefficient'],'Interpreter','latex','FontSize', 30)
leg.FontSize = 20;
set(findall(gca, 'Type', 'Line'),'LineWidth',2); SUMMARY and demonstration of database manager

check what is available in the database :

sfs = SF_Status

% demo how to load direcly the last field computed at chapter 2:

% demo how to load direcly the impedances computed at chapter 3:

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

%
fclose(fopen('SCRIPT_DEMO_ACOUSTIQUE.success','w+'));
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER GENERATED WITH SF_Status
#################################################################################

.... CONTENT OF DIRECTORY ./WORK/MESHES :
(list of meshes previously created/adapted ; couples of .msh/.ff2m files )

Index | Name              | generation mode | Date                 | Nv (nb of vertices)
1     | FFMESH_000001.msh | initial         | 09-nov.-2020 20:32:05 | 1359
2     | FFMESH_000002.msh | adapted         | 09-nov.-2020 20:32:07 | 3128
3     | FFMESH_000003.msh | initial         | 09-nov.-2020 20:33:33 | 1359
4     | FFMESH_000004.msh | adapted         | 09-nov.-2020 20:33:34 | 1061
5     | FFMESH_000005.msh | adapted         | 09-nov.-2020 20:33:37 | 2234

REMINDER : PROCEDURE TO LOAD A PREVIOUSLY COMPUTED MESH
simply type the following command  (where index is the number of the desired mesh or keyword 'last')

#################################################################################
#################################################################################

.... CONTENT OF DIRECTORY ./WORK/FORCEDFLOWS
(couples of .txt/.ff2m files )

Index | Name              | Type         | Date                 | Mesh file         | Lambda
1     | FFDATA_000001.txt | ForcedFlow   | 09-nov.-2020 20:32:06 | FFMESH_000001.msh | 2
2     | FFDATA_000002.txt | ForcedFlow   | 09-nov.-2020 20:33:34 | FFMESH_000003.msh | 1
3     | FFDATA_000003.txt | ForcedFlow   | 09-nov.-2020 20:33:36 | FFMESH_000004.msh | 0.5
4     | FFDATA_000004.txt | ForcedFlow   | 09-nov.-2020 20:33:36 | FFMESH_000004.msh | 2
5     | FFDATA_000005.txt | ForcedFlow   | 09-nov.-2020 20:33:40 | FFMESH_000005.msh | 1

REMINDER : PROCEDURE TO RECOVER A FIELD FROM THIS DIRECTORY
simply type the following command  (where index is the number of the desired field or keyword 'last')

#################################################################################

.... CONTENT OF DIRECTORY ./WORK/IMPEDANCES
(couples of .txt/.ff2m files )

Index | Name              | Type         | Date                 | Mesh file
1     | FFDATA_000001.txt | ForcedLinear | 09-nov.-2020 20:35:26 | FFMESH_000005.msh

REMINDER : PROCEDURE TO RECOVER A FIELD FROM THIS DIRECTORY
simply type the following command  (where index is the number of the desired field or keyword 'last')

#################################################################################

sfs =

struct with fields:

MESHES: [1x5 struct]
FORCEDFLOWS: [1x5 struct]
IMPEDANCES: [1x1 struct]

Forced =

struct with fields:

filename: './WORK/FORCEDFLOWS/FFDATA_000005.txt'
DataDescription: 'FORCED LINEAR RESPONSE for an axisymmetric acoustic '
datatype: 'ForcedFlow'
datastoragemode: 'CxP2'
meshfilename: './WORK/MESHES/FFMESH_000005.msh'
mesh: [1x1 struct]
Lambda: 1
INDEXING: [1x1 struct]
uz: [2234x1 double]
ur: [2234x1 double]
p: [8434x1 double]
Z: 0.4231 - 0.4727i

Imp =

struct with fields:

filename: './WORK/IMPEDANCES/FFDATA_000001.txt'
DataDescription: 'Impedance of a axisymmetric acoustic flow'
datatype: 'ForcedLinear'
datastoragemode: 'columns'