% COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $) flclear fem % COMSOL version clear vrsn vrsn.name = 'COMSOL 3.5'; vrsn.ext = ''; vrsn.major = 0; vrsn.build = 494; vrsn.rcs = '$Name: $'; vrsn.date = '$Date: 2008/09/19 16:09:48 $'; fem.version = vrsn; % Geometry g1=rect2(0.24,0.24,'base','center','pos',[-0.13,0.07]); g2=rect2('5.3/1000','5.3/1000','base','center','pos',{'0','0'},'rot','0'); g3=ellip2(5.5E-4,3.0E-4,'base','corner','pos',[-8.0E-4,0]); g4=ellip2('0.8/1000','0.8/1000','base','center','pos',{'0','0'},'rot','0'); g5=geomcomp({g2,g4},'ns',{'R1','E1'},'sf','R1-E1','edge','none'); % Analyzed geometry clear s s.objs={g5}; s.name={'CO1'}; s.tags={'g5'}; fem.draw=struct('s',s); fem.geom=geomcsg(fem); % Constants fem.const = {'R','5.3/1000', ... 'a','0.8/1000'}; % COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $) %%% Loop for the wave vector k n=10 %Number of k values for ii=1:n ii k(ii)=(ii-1)/(n-1)*0.5 % Constants fem.const = {'d','5.3/1000', ... 'r','0.8/1000', ... 'k',num2str(k(ii)), ... 'k1','1', ... 'k2','0', ... 'a1x','d', ... 'a1y','0', ... 'a2x','0', ... 'a2y','d'}; % COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $) % Geometry g5=move(g5,[5.3/1000/2,5.3/1000/2]); % Analyzed geometry clear s s.objs={g5}; s.name={'CO1'}; s.tags={'g5'}; fem.draw=struct('s',s); fem.geom=geomcsg(fem); % COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5 (COMSOL 3.5.0.494, $Date: 2008/09/19 16:09:48 $) % Initialize mesh fem.mesh=meshinit(fem, ... 'hauto',5); % Refine mesh fem.mesh=meshrefine(fem, ... 'mcase',0, ... 'rmethod','regular'); % Refine mesh fem.mesh=meshrefine(fem, ... 'mcase',0, ... 'rmethod','regular'); % Refine mesh fem.mesh=meshrefine(fem, ... 'mcase',0, ... 'rmethod','regular'); % (Default values are not included) % Application mode 1 clear appl appl.mode.class = 'AcoPressure'; appl.module = 'ACO'; appl.assignsuffix = '_acpr'; clear prop prop.analysis='eigen'; appl.prop = prop; clear bnd bnd.type = 'SH'; bnd.symmetryx = {{}}; bnd.symmetryy = {{}}; bnd.fartype = {{}}; bnd.symtypey = {{}}; bnd.symtypex = {{}}; bnd.ind = [1,1,1,1,1,1,1,1]; appl.bnd = bnd; clear equ equ.rho = 1000; equ.cs = 1500; equ.ind = [1]; appl.equ = equ; fem.appl{1} = appl; fem.frame = {'ref'}; fem.border = 1; clear units; units.basesystem = 'SI'; fem.units = units; % Scalar expressions fem.expr = {'kx','k*(k1*b1x+k2*b2x)', ... 'ky','k*(k1*b1y+k2*b2y)', ... 'b1x','2*pi*a2y/(a1x*a2y-a1y*a2x)', ... 'b1y','-2*pi*a2x/(a1x*a2y-a1y*a2x)', ... 'b2x','-2*pi*a1y/(a1x*a2y-a1y*a2x)', ... 'b2y','2*pi*a1x/(a1x*a2y-a1y*a2x)'}; % Descriptions clear descr descr.expr= {'kx','1st k-component for periodic condition','ky','2nd k-component for periodic condition','b1y','y-component of 1st reciprocal lattice vector','b2x','x-component of 2nd reciprocal lattice vector','b2y','y-component of 2nd reciprocal lattice vector','b1x','x-component of 1st reciprocal lattice vector'}; fem.descr = descr; % Coupling variable elements clear elemcpl % Extrusion coupling variables clear elem elem.elem = 'elcplextr'; elem.g = {'1'}; src = cell(1,1); clear bnd bnd.expr = {{'p',{},{}},{{},'p',{}}}; bnd.map = {{'1','1','1'},{'1','1','1'}}; bnd.ind = {{'1'},{'2'},{'3','4','5','6','7','8'}}; src{1} = {{},bnd,{}}; elem.src = src; geomdim = cell(1,1); clear bnd bnd.map = {{{},{},'2'},{{},'3',{}}}; bnd.ind = {{'1','2','5','6','7','8'},{'3'},{'4'}}; geomdim{1} = {{},bnd,{}}; elem.geomdim = geomdim; elem.var = {'pconstr1','pconstr2'}; map = cell(1,3); clear submap submap.type = 'unit'; map{1} = submap; clear submap submap.type = 'linear'; submap.sg = '1'; submap.sv = {'7','8'}; submap.dg = '1'; submap.dv = {'1','2'}; map{2} = submap; clear submap submap.type = 'linear'; submap.sg = '1'; submap.sv = {'2','8'}; submap.dg = '1'; submap.dv = {'1','7'}; map{3} = submap; elem.map = map; elemcpl{1} = elem; % Point constraint variables (used for periodic conditions) clear elem elem.elem = 'elpconstr'; elem.g = {'1'}; clear bnd bnd.constr = {{'0','pconstr2-(p*exp(j*ky*y))'},{'pconstr1-(p*exp(j*kx*x))', ... '0'}}; bnd.cpoints = {{'2','2'},{'2','2'}}; bnd.ind = {{'3'},{'4'}}; elem.geomdim = {{{},bnd,{}}}; elemcpl{2} = elem; fem.elemcpl = elemcpl; % Descriptions descr = fem.descr; descr.const= {'a2y','y component of 2st lattice vector','d','periodicidad','a1y','y component of 1st lattice vector','a1x','x component of 1st lattice vector','r','radio','k','Fraction of wave vector magnitude','k1','1st component of wave direction vector','a2x','x component of 2nd lattice vector','k2','2nd component of wave direction vector'}; fem.descr = descr; % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem); % Extend mesh fem.xmesh=meshextend(fem); % Solve problem fem.sol=femeig(fem, ... 'solcomp',{'p'}, ... 'outcomp',{'p'}, ... 'blocksize','auto', ... 'eigref','-i*2*pi*0000'); % Save current fem structure for restart purposes fem0=fem; %%% Save eigenfrequencies freky(1:6,ii)=real(fem.sol.lambda(1:6)*i/2/pi) % Plot solution postplot(fem, ... 'tridata',{'abs(p)','cont','internal','unit','Pa'}, ... 'trimap','jet(1024)', ... 'solnum',1, ... 'title','freq_acpr(1)=-5.695122e-6+21.190758i Surface: abs(p) [Pa]');%, ... %%% 'axis',[-0.0017892932446615213,0.007044900157276357,-1.4675477294442702E-4,0.005824094536405545]); end %%% Save eigenfrequencies d=str2num(cell2mat(fem.const(2))); %%% Plot dispersion relation plot(max(freky,[],1))