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:
Where is the outlet impedance which is approximated as follows :
(Ref. Fletcher & Rossing)
This approximation is equivalent to considering a reflection coefficient at the outlet as follows: (in pressure-amplitude) as follows: where
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
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 : 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 ; 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+'));