%
% 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