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