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 $\phi$ such as :

Variational formulation :

$$ \int \int_\Omega \left( \nabla \phi \cdot \nabla \phi^* + k^2 \phi \phi^*\right) dV
+ \int_{\Gamma_{in}} \phi^* dS
+ \int_{\Gamma_{out}} (i k +1/R) \phi \phi^* dV
= 0 $$

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

Contents

initialisation

close all
addpath([fileparts(fileparts(pwd)),'/SOURCES_MATLAB']);
SF_Start('verbosity',2);

set(groot, 'defaultAxesTickLabelInterpreter','latex');  % options for figures
set(groot, 'defaultLegendInterpreter','latex'); %idem
set(0,'defaultAxesFontSize',20); % idem
mkdir('html');
NOTICE  - Initialization already done
Warning: Directory already exists. 

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');
ffmesh = SF_Adapt(ffmeshInit,Forced,'Hmax',2); % Adaptation du maillage
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');
ffmesh = SF_Adapt(ffmesh,Forced1,Forced2,'Hmax',2); % Adaptation du maillage
%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 $Z(k)$

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'
    datadescriptors: 'ind,omega_r,omega_i,Z_r,Z_i,R'
       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 $Z(k)$

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:
Forced = SF_Load('FORCEDFLOWS','last')

% demo how to load direcly the impedances computed at chapter 3:
Imp = SF_Load('IMPEDANCES',1)

% [[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')
 
    mesh = SF_Load('MESH',index)
 
#################################################################################
#################################################################################
 
.... 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')
 
    object = SF_Load('FORCEDFLOWS',index)
 
#################################################################################
 
.... 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')
 
    object = SF_Load('IMPEDANCES',index)
 
#################################################################################

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'
    datadescriptors: 'p'
       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'
    datadescriptors: 'ind,omega_r,omega_i,Z_r,Z_i,R'
       meshfilename: './WORK/MESHES/FFMESH_000005.msh'
               mesh: [1x1 struct]
                ind: [401x1 double]
              omega: [401x1 double]
                  Z: [401x1 double]
                  R: [401x1 double]