Acoustic field in a pipe with harmonic forcing at the bottom

This script solves the same problem as in SCRIPT_PIPE.m and compares with theoretical predictions for the frequency.

The theoretical prediction is as follows:

$$Z_{IN} = Z_0 \frac{Z_L \cos k L + i Z_0 \sin k L}{i Z_L \sin k L +  Z_0 \cos k L}$$

Where $Z_L$ is the outlet impedance which is approximated as follows : $$ Z_L \approx Z_0 ( \frac{k^2 R^2}{2} + \frac{ 8 i k R}{3 \pi} $$

(Ref. Fletcher & Rossing)

This approximation is equivalent to considering a reflection coefficient at the outlet as follows: (in pressure-amplitude) as follows: $$ R \approx - e^{-k^2 R^2} e^{2 i k \Delta} $$ where $\Delta = 8 R/3 \pi \approx 0.85 R$

Contents

Chapter 0 : initialisation

addpath([fileparts(fileparts(pwd)) '/SOURCES_MATLAB/']);
SF_Start('verbosity',4);
NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

Chapter 1 : building an adapted mesh

ffmeshInit = SF_Mesh('Mesh_1.edp','Params',10,'problemtype','acousticaxi');
Forced = SF_LinearForced(ffmeshInit,'omega',2);
ffmesh = SF_Adapt(ffmeshInit,Forced,'Hmax',2); % Adaptation du maillage
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'  
ERROR   - SF_core_freefem: Error while using FreeFem++  ; file = Mesh_1.edp ; error code = 1 : Syntax error in your .edp program
WARNING - An error is encountered while publishing script. Launching SF_Status for diagnostics
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    ./WORK/
#################################################################################
 
.... NO MESH FOUND IN DIRECTORY ./WORK/MESHES :
 
 
#################################################################################
 
#################################################################################
Error using SF_core_log
SF_core_freefem: Error while using FreeFem++  ; file = Mesh_1.edp ; error code
= 1 : Syntax error in your .edp program

Error in SF_core_freefem (line 420)
     SF_core_log(errortype, ['SF_core_freefem: Error while using FreeFem++  ; file = ', cmd,' ; error code = ',num2str(s) ,' : ' ,errormsg]);

Error in SF_Mesh (line 67)
     value = SF_core_freefem(meshfile,'parameters',stringparam,'arguments',ffargument);

Error in SCRIPT_PIPE_Comparison_theory (line 28)
ffmeshInit = SF_Mesh('Mesh_1.edp','Params',10,'problemtype','acousticaxi');

Chapter 3 : loop over k to compute the impedance $Z(k)$

omegarange = [0.01:.01:3];

IMP = SF_LinearForced(ffmesh,omegarange)

Chapter 3b : compute theoretical results

om1 = IMP.omega;
Delta = 8/(3*pi);
L = 10;
R = exp(2i*Delta*om1).*exp(-om1.^2);
Ztheo = 1/pi*((1-R.*exp(L*2i*om1))./(1+R.*exp(L*2i*om1)));
Ztheo0 = 1/pi*((1-exp(L*2i*om1))./(1+exp(L*2i*om1))); % neglecting radiation

Plot $Z(k)$ : comparison with asymptotic

figure;
plot(IMP.omega,real(IMP.Z),'r',IMP.omega,imag(IMP.Z),'b');
hold on;
plot(om1,real(Ztheo),'r--',om1,imag(Ztheo),'b--');
legend('Z_r (num)','Z_i (num)','Z_r (asympt)', 'Z_i (asympt)')
%plot(om1,real(Ztheo0),'r:',om1,imag(Ztheo0),'b:');

title(['Impedance $Z_r$ and $Z_i$, comparison with approx.'],'Interpreter','latex','FontSize', 24)
xlabel('$\omega R/c$','Interpreter','latex','FontSize', 24);
ylabel('$Z_r,Z_i$','Interpreter','latex','FontSize', 24);
set(findall(gca, 'Type', 'Line'),'LineWidth',2);
pause(0.1);

figure;
semilogy(IMP.omega,abs(IMP.Z),'r');
hold on;
semilogy(om1,abs(Ztheo),'r:');
title(['Impedance $|Z|$ , comparision with asympt. approx'],'Interpreter','latex','FontSize', 24)
xlabel('$\omega R/c$','Interpreter','latex','FontSize', 24);
ylabel('$|Z|$','Interpreter','latex','FontSize', 24);
legend('|Z| (num)', '|Z| (theory)')
set(findall(gca, 'Type', 'Line'),'LineWidth',2);
pause(0.1);

Plot $Z(k)$ ; looking for asymptotes

figure;
plot(IMP.omega,real(IMP.Z)./(IMP.omega).^2,'r',IMP.omega,imag(IMP.Z)./IMP.omega,'b--');
xlim([0  .08]);
title(['Low-freq. asymptotic behaviours for $Z_r$ and $Z_i$'],'Interpreter','latex','FontSize', 24)
xlabel('$\omega R/c$','Interpreter','latex','FontSize', 24);
ylabel('$Z_r/\omega^2 ,Z_i/\omega$','Interpreter','latex','FontSize', 24);
om0 = [0:.001:.08];
hold on; plot(om0,1/(2*pi)*ones(size(om0)),'r:',om0,-(10+8/(3*pi))/pi*ones(size(om0)),'b:');
set(findall(gca, 'Type', 'Line'),'LineWidth',2);
legend('Z_r (num)','Z_i (num)','Z_r (asympt)', 'Z_i (asympt)')
pause(0.1);

plot Reflection coefficient

figure;
plot(IMP.omega,IMP.R,'b-');
xlabel('$\omega R/c$','Interpreter','latex','FontSize', 24);
ylabel('$R$','Interpreter','latex','FontSize', 24);
title(['Reflection coefficient (energy)'],'Interpreter','latex','FontSize', 24)
hold on; plot(IMP.omega,exp(-2*IMP.omega.^2),'b--');
legend('num','theory')

set(findall(gca, 'Type', 'Line'),'LineWidth',2);

[[PUBLISH]] this tag is to allow automatic publication as html

%
fclose(fopen('SCRIPT_PIPE_Comparison_theory.success','w+'));