### FEM Simulation of the Toneholes Studied by Dalmont et al. (2002)

toneholeflangeddalmont.m
```%
% COMSOL / Matlab script
% FEM Simulation of two toneholes with the same geometry as in
% Dalmont et al. (2002), Experimental determination of the equivalent circuit
% of an open side hole: Linear and non linear behavior. Acustica, 88, 567-575.
%
% Antoine Lefebvre, McGill University
%
% Copyright (c) 2010 Antoine Lefebvre
%
flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 603;
vrsn.rcs = '\$Name:  \$';
vrsn.date = '\$Date: 2008/12/03 17:02:19 \$';
fem.version = vrsn;

rho = 1.184;
c   = 346.16; % speed of sound
a   = 0.01; % radius of the body
S   = pi*a^2;
Zc  = rho*c/S;

ka = [0.01 0.1:0.05:1.00];

k = ka/a;
f = k*c/(2*pi);

L = 0.05; % half length of the body

cases{1} = [0.007, 0.0092]; %b, t
cases{2} = [0.010, 0.0101]; %b, t

R = 0.15; % radius of the radiation sphere

for i=1:2
b = cases{i}(1);
t = cases{i}(2);

mesh_edge_size = 2*pi*b/100;

S_hole = pi*b^2;
Zc_hole = rho*c/S_hole;

filename = sprintf('case%d.mat',i);

% Geometry
body=cylinder3(a,2*L,'pos',[-L,0,0],'axis',[1,0,0],'rot','0');
tonehole=cylinder3(b,a+t);
instrument=geomdel(body+tonehole);
exterior=sphere3(R,'pos',[0,0,a+t]);
block=block3(2*R,2*R,R,'base','corner','pos',[-R,-R,a+t]);
model=instrument+exterior*block;
halfmodel=model*block3(2*R,2*R,2*a+t+R,'base','corner','pos',[-R,0,-a]);
fourthmodel=halfmodel*block3(R,R,2*a+t+R,'base','corner','pos',[-R,0,-a]);

% Analyzed geometry
clear s
s.objs={fourthmodel};
s.name={'CO1'};
s.tags={'halfmodel'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

fem.const = {'rho_air','1.184 [kg/m^3]', 'cs_air','346.16 [m/s]', ...
'Zc','rho_air*cs_air/S', 'a','0.01 [m]', 'S','pi*a^2', ...
'L','0.05 [m]'};

fem.mesh=meshinit(fem, 'hauto',5, 'hmaxedg', ...
[11,mesh_edge_size,13,mesh_edge_size]);

clear appl
appl.mode.class = 'AcoPressure';
appl.module = 'ACO';
appl.assignsuffix = '_acpr';

clear prop
prop.elemdefault='Lag3';
appl.prop = prop;

%%%% symmetric boundary condition (normal acceleration is null)
clear bnd_s
bnd_s.p0 = {0,0,0,1,0};
bnd_s.type = {'SH','RAD','NA','RAD','cont'};
bnd_s.wavetype = {'PL','SPH','PL','PL','PL'};
bnd_s.ind = [3,1,2,4,3,1,1,1,5,3,3];

appl.bnd = bnd_s;

%%%% anti-symmetric boundary condition (pressure is null)
clear bnd_a
bnd_a.p0 = {0,0,0,1,0,0};
bnd_a.type = {'SH','RAD','NA','RAD','cont','SS'};
bnd_a.wavetype = {'PL','SPH','PL','PL','PL','PL'};
bnd_a.ind = [3,1,2,4,3,1,1,1,5,6,6];

clear equ
equ.rho = 'rho_air';
equ.cs = 'cs_air';
equ.ind = [1,1];
appl.equ = equ;
appl.var = {'freq','freq'};
fem.appl{1} = appl;
fem.border = 1;

fem=multiphysics(fem);
fem.xmesh=meshextend(fem);

fem.sol=femstatic(fem, 'solcomp',{'p'}, 'outcomp',{'p'}, 'pname','freq', ...
'plist',f, 'oldcomp',{},'linsolver','gmres');

pin_s   =  postint(fem,'2*p/S',                     'Dl',[4],'Edim',[2],...
'Solnum',[1:length(f)]);
ZcUin_s = -postint(fem,'2*rho_air*cs_air*nv_acpr/S','Dl',[4],'Edim',[2],...
'Solnum',[1:length(f)]);

fem.appl{1}.bnd = bnd_a;
fem=multiphysics(fem);
fem.xmesh=meshextend(fem);
fem.sol=femstatic(fem, 'solcomp',{'p'}, 'outcomp',{'p'}, 'pname','freq',...
'plist',f, 'oldcomp',{},'linsolver','gmres');

pin_a   =  postint(fem,'2*p/S',                     'Dl',[4],'Edim',[2],...
'Solnum',[1:length(f)]);
ZcUin_a = -postint(fem,'2*rho_air*cs_air*nv_acpr/S','Dl',[4],'Edim',[2],...
'Solnum',[1:length(f)]);

TH = CalculateTMfromSymmetricResults(k*L, pin_s, ZcUin_s, pin_a, ZcUin_a);

Zs = Zc./TH(3,:);
Za = 2*( TH(1,:)-1) .* Zs;

ta = real(Za./(j*k*Zc_hole));
ts = real(Zs./(j*k*Zc_hole));
shunt_losses = real(Zs./Zc_hole)./((k*b).^2);

delta=b/a;
save(filename,'t','b','a','delta','k','Zs','Za','ts','ta','TH','Zc',...
'Zc_hole','shunt_losses');
end```