% % COMSOL / Matlab script % Calculation of the acoustic radiation impedance of an unflanged pipe % with the FEM (Helmholtz equation) % % Antoine Lefebvre, McGill University % % Copyright (c) 2010 Antoine Lefebvre % flclear fem rho = 1.25; c = 343; a = 0.01; L = 0.03; N = 5; ka = [0.01 colon(0.1,0.1,1.0)]; k = ka/a; block = block3(2*N*a,2*N*a,2*N*a, 'pos',[0,0,-N*a]); sphere = sphere3(N*a)*block; pipe = cylinder3(a,L,'pos', [0,0,-L])*block; radiation = sphere-pipe; % Analyzed geometry clear s s.objs={radiation,pipe}; s.name={'rad','pipe'}; s.tags={'rad','pipe'}; fem.draw=struct('s',s); fem.geom=geomgroup(fem,'imprint','off','paircand','none'); edge=0.001; fem.mesh=meshinit(fem, 'hmax',0.005, 'hmaxfact',0.8, 'hcurve',0.5, ... 'hgrad',1.3, 'hcutoff',0.02, 'hnarrow',0.6, ... 'hmaxedg', [13,edge,23,edge], 'point',[], ... 'edge',[], 'face',[], 'subdomain',[2]); fem.mesh=meshcopy(fem, 'source',11, 'target',5, 'mcase',0); fem.mesh=meshinit(fem, 'hmax',0.025, 'hmaxfact',0.8, 'hcurve',0.5, ... 'hgrad',1.3, 'hcutoff',0.02, 'hnarrow',0.6, ... 'hmaxedg', [13,edge,23,edge], 'point',[], ... 'edge',[], 'face',[], 'subdomain',[1], ... 'meshstart',fem.mesh); % Constants fem.const = {'v0','1 [m/s]', 'a','0.01 [m]', 'S','pi*a^2', 'c','343 [m/s]', ... 'rho','1.25 [kg/m^3]'}; clear appl appl.mode.class = 'AcoPressure'; appl.module = 'ACO'; appl.sshape = 3; appl.assignsuffix = '_acpr'; clear prop prop.elemdefault='Lag3'; appl.prop=prop; clear bnd bnd.p0 = {0,0,1,0}; bnd.type = {'SH','RAD','RAD','NA'}; bnd.wavetype = {'PL','SPH','PL','PL'}; bnd.ind = [4,4,2,1,1,2,1,4,4,3,1,1]; clear pair pair.type = 'cont'; pair.pair = 'open end'; bnd.pair = pair; appl.bnd = bnd; appl.var = {'freq','freq'}; fem.appl{1} = appl; fem.border = 1; % Boundary pairs clear pair pair{1}.type= 'identity'; pair{1}.name= 'open end'; pair{1}.src.dl = [11]; pair{1}.src.operator = 'src2dst_ip1'; pair{1}.dst.dl = [5]; pair{1}.dst.operator = 'dst2src_ip1'; bnd.pair = pair; fem.bnd = bnd; fem.expr = {'freq','ka*c/(2*pi*a)','a0','i*omega_acpr*v0'}; fem=multiphysics(fem); fem.xmesh=meshextend(fem); fem.sol=femstatic(fem, 'solcomp',{'p'}, 'outcomp',{'p'}, 'blocksize','auto', ... 'pname','ka', 'plist',ka, 'oldcomp',{}); pin = postint(fem,'p/S','Dl',[10],'Edim',[2],'Solnum',[1:length(ka)]); ZcUin = -postint(fem,'rho*c*nv_acpr/S','Dl',[10],'Edim',[2],'Solnum',[1:length(ka)]); Zin = pin./ZcUin; Zr = (Zin.*cos(k*L)-i*sin(k*L))./(-i*Zin.*sin(k*L)+cos(k*L)); flsave('unflanged_pipe_quarter.mph', fem); save('unflanged_pipe_quarter.mat', 'Zr','a','ka');