Note: This discussion is about an older version of the COMSOL Multiphysics® software. The information provided may be out of date.
Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.
I get different results when using different mesh parameters 4 or 5 using script file below for modelling plasmonic structure
Posted 20 janv. 2011, 22:17 UTC−5 RF & Microwave Engineering, Mesh, Studies & Solvers Version 3.5a 0 Replies
Please login with a confirmed email address before reporting spam
Can anyone help with the problem of getting different transmission curves, when modelling plasmon resonances using Comsol 3.5 and a Matlab script as below, when applying Mesh parameter 4 and mesh parameter 5. How to avoid problems like that.
The solution with Mesh parameter 4 and a Mesh 5 has a common plasmon peak around 830 nm, but differes in other regions.
Please help,
Krystyna.
%geometry definition:
R=35; %Cylinder radius in nm
H=40; %Cylinder height in nm
r=5; %Cylinder corner radius in nm
L=150; %Unit cell side in nm
H1=50; %Glass prism height in nm
H2=100; %Air height in nm
%H3=15; %Protective layer height in nm
%H4=30; %mirror height in nm
%Accuracy:
ACC=5; % Between 1 and 9, the lower the finer mesh
%Calculations:
flclear fem
a=rect2(R,H);
a=fillet(a,'point',[2 3],'radii',r);
a=revolve(a);
a=rotate(a,pi/2,[1 0 0]);
g2=geomcopy({a});
g3=move(g2,[-35,0,0]);
g4=geomcopy({g3});
g5=move(g4,[70,0,0]);
b=block3(L,L,H1,'pos',[-L/2 -L/2 -H1]);
%b=block3(L,L,H1,'pos',[-L/2 -L/2 -H1-H3]);
c=block3(L,L,H2,'pos',[-L/2 -L/2 0]);
%d=block3(L,L,H3,'pos',[-L/2 -L/2 -H3]);
%e=block3(L,L,H4,'pos',[-L/2 -L/2 -H3-H4]);
fem.geom=scale(g3+g5+b+c,1e-9,1e-9,1e-9);
%fem.geom=scale(a+b+c+d+e,1e-9,1e-9,1e-9);
%geomplot(fem,'facelabels','on','transparency',0.1)
%xlabel('x')
%ylabel('y')
%zlabel('z')
%return
%mesh creation:
%hauto = the size of the mesh, 1=finest, 9=coarsest
fem.mesh=meshinit(fem,'face',[1 4 ],'subdomain','none','hauto',ACC);
fem.mesh=meshcopy(fem,'source',1,'target',50);
fem.mesh=meshcopy(fem,'source',4,'target',51);
fem.mesh=meshinit(fem,'face',[2 5],'subdomain','none','hauto',ACC,'meshstart',fem.mesh);
fem.mesh=meshcopy(fem,'source',2,'target',8);
fem.mesh=meshcopy(fem,'source',5,'target',9);
%meshplot(fem,'boundmode','on','transparency',0.5)
%return
fem.mesh=meshinit(fem,'hauto',ACC,'meshstart',fem.mesh);
%meshplot(fem,'sdl',5)
%return
%filling the application with structure:
fem.appl{1}.mode.class='ElectromagneticWaves';
fem.appl{1}.module='RF';
fem.appl{1}.gporder=4;
fem.appl{1}.cporder=2;
fem.appl{1}.assignsuffix='_rfw';
fem.appl{1}.prop.inputvar='lambda';
fem.appl{1}.bnd.type={'SC' 'cont' 'periodic' 'periodic' 'SC'};
fem.appl{1}.bnd.E0={{0;0;0} {0;0;0} {0;0;0} {0;0;0} {'cos(theta)';0;'sin(theta)'}};
fem.appl{1}.bnd.index={[0] [0] [1] [0] [0]};
fem.appl{1}.bnd.kper={{0;0;0} {0;0;0} {'k0*sin(theta)';0;'-k0*cos(theta)'} {0;0;0} {0;0;0}};
fem.appl{1}.bnd.pertype={'sym' 'sym' 'floque' 'sym' 'sym'};
%fem.appl{1}.bnd.ind=[3 4 1 3 4 2 5 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3];
fem.appl{1}.bnd.ind=[3,4,1,3,4,2,5,4,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, ...
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3];
% 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
% 0 1 2 3 4
fem.appl{1}.equ.n={[1.55] [1] 'nAgr(lambda)-i*nAgi(lambda)'};
fem.appl{1}.equ.matparams='n';
fem.appl{1}.equ.ind=[1 2 3 3];
% 1 2 3 4 5
fem.appl{1}.var={'lambda0' 'lambda'};
%defining the functions:
fem.functions{1}.type='interp';
fem.functions{1}.name='nAgr';
fem.functions{1}.method='piecewisecubic';
fem.functions{1}.extmethod='interior';
fem.functions{1}.x={'1.8785e-007' '1.9163e-007' '1.9525e-007' '1.9933e-007' '2.0325e-007' '2.0733e-007' '2.1194e-007' '2.1638e-007' '2.2140e-007' '2.2625e-007' '2.3131e-007' '2.3706e-007' '2.4263e-007' '2.4896e-007' '2.5511e-007' '2.6157e-007' '2.6895e-007' '2.7613e-007' '2.8437e-007' '2.9242e-007' '3.0093e-007' '3.1074e-007' '3.2037e-007' '3.3151e-007' '3.4250e-007' '3.5424e-007' '3.6791e-007' '3.8149e-007' '3.9739e-007' '4.1328e-007' '4.3351e-007' '4.5085e-007' '4.7142e-007' '4.9594e-007' '5.2094e-007' '5.4860e-007' '5.8209e-007' '6.1684e-007' '6.5949e-007' '7.0446e-007' '7.5600e-007' '8.2109e-007' '8.9197e-007' '9.8400e-007' '1.0876e-006' '1.2155e-006' '1.3931e-006' '1.6102e-006' '1.9373e-006'};
fem.functions{1}.data={'1.07' '1.10' '1.12' '1.14' '1.15' '1.18' '1.20' '1.22' '1.25' '1.26' '1.28' '1.28' '1.30' '1.31' '1.33' '1.35' '1.38' '1.41' '1.41' '1.39' '1.34' '1.13' '0.61' '0.17' '0.14' '0.10' '0.07' '0.05' '0.05' '0.05' '0.04' '0.04' '0.05' '0.05' '0.05' '0.06' '0.05' '0.06' '0.05' '0.04' '0.03' '0.04' '0.04' '0.04' '0.04' '0.09' '0.13' '0.15' '0.24'};
fem.functions{2}.type='interp';
fem.functions{2}.name='nAgi';
fem.functions{2}.method='piecewisecubic';
fem.functions{2}.extmethod='interior';
fem.functions{2}.x={'1.8785e-007' '1.9163e-007' '1.9525e-007' '1.9933e-007' '2.0325e-007' '2.0733e-007' '2.1194e-007' '2.1638e-007' '2.2140e-007' '2.2625e-007' '2.3131e-007' '2.3706e-007' '2.4263e-007' '2.4896e-007' '2.5511e-007' '2.6157e-007' '2.6895e-007' '2.7613e-007' '2.8437e-007' '2.9242e-007' '3.0093e-007' '3.1074e-007' '3.2037e-007' '3.3151e-007' '3.4250e-007' '3.5424e-007' '3.6791e-007' '3.8149e-007' '3.9739e-007' '4.1328e-007' '4.3351e-007' '4.5085e-007' '4.7142e-007' '4.9594e-007' '5.2094e-007' '5.4860e-007' '5.8209e-007' '6.1684e-007' '6.5949e-007' '7.0446e-007' '7.5600e-007' '8.2109e-007' '8.9197e-007' '9.8400e-007' '1.0876e-006' '1.2155e-006' '1.3931e-006' '1.6102e-006' '1.9373e-006'};
fem.functions{2}.data={'1.168' '1.203' '1.226' '1.251' '1.277' '1.304' '1.330' '1.387' '1.427' '1.460' '1.497' '1.536' '1.577' '1.631' '1.688' '1.749' '1.803' '1.847' '1.869' '1.878' '1.889' '1.893' '1.898' '1.883' '1.871' '1.866' '1.895' '1.933' '1.952' '1.958' '1.948' '1.914' '1.849' '1.833' '2.081' '2.455' '2.863' '3.272' '3.697' '4.103' '4.542' '5.083' '5.663' '6.350' '7.150' '8.145' '9.519' '11.210' '13.780'};
%defining the constants:
fem.const={'lambda' '500e-9' 'theta' 'angle*pi/180' 'angle' '0' 'k0' '2*pi/lambda'};
%creating the multiphysics application:
fem=multiphysics(fem);
%solve the problem:
fem=femdiff(fem);
fem.xmesh=meshextend(fem);
%Alternative 1, just solve for one wavelength only and plot the field:
fem.sol=femlin(fem,'linsolver','PARDISO');
figure
postslice(fem,'normE_rfw','slicexspacing',0,'sliceyspacing',1,'slicezspacing',0);
%Alternative 2, sweep over wavelength and plot the power out or absorbed power versus wavelength:
lambda=[350:10:1100]*1e-9;
fem.sol=femlin(fem,'linsolver','PARDISO','pname',{'lambda'},'plist',lambda);
%Power out:
Pout=postint(fem,'nPoav_rfw','edim',2,'dl',[3],'solnum','all')
figure
plot(lambda/1e-9,Pout/1e-15,'k','linewidth',2);
set(gca,'linewidth',2,'fontname','times new roman','fontsize',12);
set(gcf,'color','w');
xlabel('\lambda (nm)');
ylabel('Power out (fW)');
%Absorbed power:
Pabs= postint(fem,'Qav_rfw','edim',3,'dl',[3 4],'solnum','all')
figure
plot(lambda/1e-9,-Pres/1e-15,'k','linewidth',2);
set(gca,'linewidth',2,'fontname','times new roman','fontsize',12);
set(gcf,'color','w');
xlabel('\lambda (nm)');
ylabel('Absorbed power (fW)');
set(gca,'xlim',lambda([1 end])/1e-9)
The solution with Mesh parameter 4 and a Mesh 5 has a common plasmon peak around 830 nm, but differes in other regions.
Please help,
Krystyna.
%geometry definition:
R=35; %Cylinder radius in nm
H=40; %Cylinder height in nm
r=5; %Cylinder corner radius in nm
L=150; %Unit cell side in nm
H1=50; %Glass prism height in nm
H2=100; %Air height in nm
%H3=15; %Protective layer height in nm
%H4=30; %mirror height in nm
%Accuracy:
ACC=5; % Between 1 and 9, the lower the finer mesh
%Calculations:
flclear fem
a=rect2(R,H);
a=fillet(a,'point',[2 3],'radii',r);
a=revolve(a);
a=rotate(a,pi/2,[1 0 0]);
g2=geomcopy({a});
g3=move(g2,[-35,0,0]);
g4=geomcopy({g3});
g5=move(g4,[70,0,0]);
b=block3(L,L,H1,'pos',[-L/2 -L/2 -H1]);
%b=block3(L,L,H1,'pos',[-L/2 -L/2 -H1-H3]);
c=block3(L,L,H2,'pos',[-L/2 -L/2 0]);
%d=block3(L,L,H3,'pos',[-L/2 -L/2 -H3]);
%e=block3(L,L,H4,'pos',[-L/2 -L/2 -H3-H4]);
fem.geom=scale(g3+g5+b+c,1e-9,1e-9,1e-9);
%fem.geom=scale(a+b+c+d+e,1e-9,1e-9,1e-9);
%geomplot(fem,'facelabels','on','transparency',0.1)
%xlabel('x')
%ylabel('y')
%zlabel('z')
%return
%mesh creation:
%hauto = the size of the mesh, 1=finest, 9=coarsest
fem.mesh=meshinit(fem,'face',[1 4 ],'subdomain','none','hauto',ACC);
fem.mesh=meshcopy(fem,'source',1,'target',50);
fem.mesh=meshcopy(fem,'source',4,'target',51);
fem.mesh=meshinit(fem,'face',[2 5],'subdomain','none','hauto',ACC,'meshstart',fem.mesh);
fem.mesh=meshcopy(fem,'source',2,'target',8);
fem.mesh=meshcopy(fem,'source',5,'target',9);
%meshplot(fem,'boundmode','on','transparency',0.5)
%return
fem.mesh=meshinit(fem,'hauto',ACC,'meshstart',fem.mesh);
%meshplot(fem,'sdl',5)
%return
%filling the application with structure:
fem.appl{1}.mode.class='ElectromagneticWaves';
fem.appl{1}.module='RF';
fem.appl{1}.gporder=4;
fem.appl{1}.cporder=2;
fem.appl{1}.assignsuffix='_rfw';
fem.appl{1}.prop.inputvar='lambda';
fem.appl{1}.bnd.type={'SC' 'cont' 'periodic' 'periodic' 'SC'};
fem.appl{1}.bnd.E0={{0;0;0} {0;0;0} {0;0;0} {0;0;0} {'cos(theta)';0;'sin(theta)'}};
fem.appl{1}.bnd.index={[0] [0] [1] [0] [0]};
fem.appl{1}.bnd.kper={{0;0;0} {0;0;0} {'k0*sin(theta)';0;'-k0*cos(theta)'} {0;0;0} {0;0;0}};
fem.appl{1}.bnd.pertype={'sym' 'sym' 'floque' 'sym' 'sym'};
%fem.appl{1}.bnd.ind=[3 4 1 3 4 2 5 4 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3];
fem.appl{1}.bnd.ind=[3,4,1,3,4,2,5,4,4,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, ...
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3];
% 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
% 0 1 2 3 4
fem.appl{1}.equ.n={[1.55] [1] 'nAgr(lambda)-i*nAgi(lambda)'};
fem.appl{1}.equ.matparams='n';
fem.appl{1}.equ.ind=[1 2 3 3];
% 1 2 3 4 5
fem.appl{1}.var={'lambda0' 'lambda'};
%defining the functions:
fem.functions{1}.type='interp';
fem.functions{1}.name='nAgr';
fem.functions{1}.method='piecewisecubic';
fem.functions{1}.extmethod='interior';
fem.functions{1}.x={'1.8785e-007' '1.9163e-007' '1.9525e-007' '1.9933e-007' '2.0325e-007' '2.0733e-007' '2.1194e-007' '2.1638e-007' '2.2140e-007' '2.2625e-007' '2.3131e-007' '2.3706e-007' '2.4263e-007' '2.4896e-007' '2.5511e-007' '2.6157e-007' '2.6895e-007' '2.7613e-007' '2.8437e-007' '2.9242e-007' '3.0093e-007' '3.1074e-007' '3.2037e-007' '3.3151e-007' '3.4250e-007' '3.5424e-007' '3.6791e-007' '3.8149e-007' '3.9739e-007' '4.1328e-007' '4.3351e-007' '4.5085e-007' '4.7142e-007' '4.9594e-007' '5.2094e-007' '5.4860e-007' '5.8209e-007' '6.1684e-007' '6.5949e-007' '7.0446e-007' '7.5600e-007' '8.2109e-007' '8.9197e-007' '9.8400e-007' '1.0876e-006' '1.2155e-006' '1.3931e-006' '1.6102e-006' '1.9373e-006'};
fem.functions{1}.data={'1.07' '1.10' '1.12' '1.14' '1.15' '1.18' '1.20' '1.22' '1.25' '1.26' '1.28' '1.28' '1.30' '1.31' '1.33' '1.35' '1.38' '1.41' '1.41' '1.39' '1.34' '1.13' '0.61' '0.17' '0.14' '0.10' '0.07' '0.05' '0.05' '0.05' '0.04' '0.04' '0.05' '0.05' '0.05' '0.06' '0.05' '0.06' '0.05' '0.04' '0.03' '0.04' '0.04' '0.04' '0.04' '0.09' '0.13' '0.15' '0.24'};
fem.functions{2}.type='interp';
fem.functions{2}.name='nAgi';
fem.functions{2}.method='piecewisecubic';
fem.functions{2}.extmethod='interior';
fem.functions{2}.x={'1.8785e-007' '1.9163e-007' '1.9525e-007' '1.9933e-007' '2.0325e-007' '2.0733e-007' '2.1194e-007' '2.1638e-007' '2.2140e-007' '2.2625e-007' '2.3131e-007' '2.3706e-007' '2.4263e-007' '2.4896e-007' '2.5511e-007' '2.6157e-007' '2.6895e-007' '2.7613e-007' '2.8437e-007' '2.9242e-007' '3.0093e-007' '3.1074e-007' '3.2037e-007' '3.3151e-007' '3.4250e-007' '3.5424e-007' '3.6791e-007' '3.8149e-007' '3.9739e-007' '4.1328e-007' '4.3351e-007' '4.5085e-007' '4.7142e-007' '4.9594e-007' '5.2094e-007' '5.4860e-007' '5.8209e-007' '6.1684e-007' '6.5949e-007' '7.0446e-007' '7.5600e-007' '8.2109e-007' '8.9197e-007' '9.8400e-007' '1.0876e-006' '1.2155e-006' '1.3931e-006' '1.6102e-006' '1.9373e-006'};
fem.functions{2}.data={'1.168' '1.203' '1.226' '1.251' '1.277' '1.304' '1.330' '1.387' '1.427' '1.460' '1.497' '1.536' '1.577' '1.631' '1.688' '1.749' '1.803' '1.847' '1.869' '1.878' '1.889' '1.893' '1.898' '1.883' '1.871' '1.866' '1.895' '1.933' '1.952' '1.958' '1.948' '1.914' '1.849' '1.833' '2.081' '2.455' '2.863' '3.272' '3.697' '4.103' '4.542' '5.083' '5.663' '6.350' '7.150' '8.145' '9.519' '11.210' '13.780'};
%defining the constants:
fem.const={'lambda' '500e-9' 'theta' 'angle*pi/180' 'angle' '0' 'k0' '2*pi/lambda'};
%creating the multiphysics application:
fem=multiphysics(fem);
%solve the problem:
fem=femdiff(fem);
fem.xmesh=meshextend(fem);
%Alternative 1, just solve for one wavelength only and plot the field:
fem.sol=femlin(fem,'linsolver','PARDISO');
figure
postslice(fem,'normE_rfw','slicexspacing',0,'sliceyspacing',1,'slicezspacing',0);
%Alternative 2, sweep over wavelength and plot the power out or absorbed power versus wavelength:
lambda=[350:10:1100]*1e-9;
fem.sol=femlin(fem,'linsolver','PARDISO','pname',{'lambda'},'plist',lambda);
%Power out:
Pout=postint(fem,'nPoav_rfw','edim',2,'dl',[3],'solnum','all')
figure
plot(lambda/1e-9,Pout/1e-15,'k','linewidth',2);
set(gca,'linewidth',2,'fontname','times new roman','fontsize',12);
set(gcf,'color','w');
xlabel('\lambda (nm)');
ylabel('Power out (fW)');
%Absorbed power:
Pabs= postint(fem,'Qav_rfw','edim',3,'dl',[3 4],'solnum','all')
figure
plot(lambda/1e-9,-Pres/1e-15,'k','linewidth',2);
set(gca,'linewidth',2,'fontname','times new roman','fontsize',12);
set(gcf,'color','w');
xlabel('\lambda (nm)');
ylabel('Absorbed power (fW)');
set(gca,'xlim',lambda([1 end])/1e-9)
0 Replies Last Post 20 janv. 2011, 22:17 UTC−5
Hello Krystyna Drozdowicz-Tomsia
Your Discussion has gone 30 days without a reply. If you still need help with COMSOL and have an on-subscription license, please visit our Support Center for help.
If you do not hold an on-subscription license, you may find an answer in another Discussion or in the Knowledge Base.