COMSOL光学应用

利用COMSOL提升研究生表面等离激元光学的教学质量作者：黄勇刚来源：《科技视界》2019年第19期【摘要】表面等離激元具有极强的局域场增强效应，能突破光学衍射极限，为纳米尺度上光的调控提供了一种可能的途径，被认为是实现纳米量级集成光子学器件的最优方式之一，近年来已发展成一门新型的学科——表面等离激元光子学。

将基于有限元的COMSOL软件引入表面等离激元光学的教学过程中，能有效地将科学前沿与研究生课程教学有机地结合在一起，用学科的前沿案例吸引学生的学习兴趣，用课程的系统性理论提供科学研究的新思路，以达到科研与教学互相促进的作用。

本文结合笔者9年的教学经验，以COMSOL在表面等离激元光波导教学中的应用为例，介绍波导模面积和传播常数等工作参量的求解过程，加深了学生对表面等离激元波导模式的理解，激发了学生学习的主动性，提高学生仿真实践能力和设计新型波导的创新能力，对提升人才培养质量具有积极作用。

【关键词】表面等离激元光波导;COMSOL;模面积;传播距离中图分类号： O485 文献标识码： A 文章编号： 2095-2457（2019）19-0067-002DOI：10.19694/ki.issn2095-2457.2019.19.0310 引言表面等离激元是金属中自由电子和电磁场共谐震荡量子化后的准粒子[1-6]。

根据自由电子振荡受到结构尺寸的限制不同，可分为局域表面等离激元（三个维度受限，如金属纳米颗粒）和传导的表面等离激元（至少一个维度不受限，如金属和介质分界面或金属纳米线等）。

表面等离激元能把光场压缩在突破传统光学衍射极限的纳米尺度，可用来减小光学器件的尺寸以及增强光与物质的相互作用等，在自发辐射增强、光催化、单分子拉曼增强、能源、生物医学高灵敏探测以及量子信息与量子计算等领域具有重要的应用，为实现全光集成，发展更小、更快和更高效的纳米光子学器件提供了一条有效途径，也为能源开发和利用效率的提高发挥作用，目前已发展成为一门新兴的学科——表面等离激元光子学[7]。

基于COMSOL的HID灯物理模型【摘要】本文介绍了基于COMSOL的HID灯物理模型。

我们阐述了研究背景和研究意义。

接着，详细解释了HID灯的工作原理以及COMSOL在光学仿真中的应用。

然后，我们描述了如何使用COMSOL建立HID灯的光学模型，并对模型进行了分析和光学参数优化。

结论部分讨论了基于COMSOL的HID灯光学模型的可靠性，对HID灯设计的启示，以及未来的展望。

这篇文章对于研究者和工程师在设计和优化HID灯方面具有重要参考价值。

【关键词】HID灯、COMSOL、光学仿真、光学模型、光学参数、优化、可靠性、设计、未来展望1. 引言1.1 研究背景研究背景：HID灯（高强度气体放电灯）是一种高效、高亮度的照明设备，广泛应用于汽车大灯、室外照明等领域。

随着LED等新光源技术的发展，HID灯的照明效果与能效比逐渐受到挑战。

对HID灯的光学性能进行深入研究和优化已成为当今照明行业的热点问题。

基于COMSOL来建立HID灯的光学模型，对于深入了解其光学性能、提高设计效率、降低成本具有重要意义。

通过模拟分析HID灯的光学参数，能够优化灯具结构，提高其光效和光质，为照明行业的发展带来新的启示。

1.2 研究意义HID灯作为目前常见的一种高效、高亮度照明设备，广泛应用于街道照明、工业照明等领域。

随着人们对照明质量和节能环保要求的不断提高，对HID灯的设计和优化也显得尤为重要。

基于COMSOL的HID灯物理模型可以帮助我们更深入地了解HID灯的光学特性，优化其光学参数，提高照明效果和能效比。

通过建立基于COMSOL的HID 灯光学模型，可以实现对HID灯光束形状、亮度分布等关键参数的精确模拟和优化，为HID灯的设计和制造提供重要参考。

研究基于COMSOL的HID灯光学模型具有重要的实际意义和应用前景，将有助于推动HID灯的发展，提高其在照明领域的竞争力。

2. 正文2.1 HID灯的工作原理HID灯(High Intensity Discharge lamp)是一种高强度放电灯，也称为气体放电灯。

此模型基于《COMSOL 软件许可协议》6.0 版本授权。

所有商标均为其各自所有者的财产。

请参见 /trademarks 。

在 COMSOL Multiphysics 6.0 版本中创建光学环形谐振腔陷波滤波器2 |光学环形谐振腔陷波滤波器简介最简单的光学环形谐振腔由直波导和环形波导组成。

这两个波导的芯层彼此靠近放置，因此光从一个波导耦合到另一个波导。

当环形波导的长度为波长的整数倍时，环形波导与波长产生谐振，并且存储在圆环中的光波能量增加。

通过直波导传输的波是入射波与从环耦合到直波导的波形成的干涉波。

您可以大致地将环形谐振腔视为如下面的图 1 所示。

入射波 E i 1 的一部分在直波导中传输，而其中的一部分场耦合到环形波导中。

类似地，环形波导中的部分光与直波导中的光相耦合，而其余的波继续围绕环形波导传输。

图 1：光学环形谐振腔示意图，显示入射场 E i1 和 E i2 以及透射场/耦合场 E t1 和 E t2，还指出了传输系数 t 和耦合系数 κ，以及往返损耗 L 。

透射场通过以下矩阵-矢量关系式与入射场相关(1)上面定义的矩阵单元确保总输入功率等于总输出功率，(2)其中假设传输系数和耦合系数之间的关系如下(3)tt *κ-κ*E i1E t1E i2E t2LE t 1E t 2t κκ*–t*E i 1E i 2=E t 12E t 22+E i 12E i 22+=t 2κ2+1=3 |光学环形谐振腔陷波滤波器此外，当波围绕环形波导传播时，您会得到以下关系式(4)其中 L 是围绕环形波导传播的损耗系数， 是累积的相位。

结合方程 1、方程 3 和方程 4，可以将透射场写为(5)这里，透射系数被分成传输损耗 |t | 和相应的相位，(6)请注意，谐振情况下，当 是 2π 的整数倍，且 |t | = L 时，透射场为零。

|t | = L 的条件称为临界耦合。

因此，当耦合器传输损耗与环形波导周围传播的波的损耗相抵消时，您会得到环形谐振腔用作带阻滤波器（即陷波滤波器）的最佳条件。

Comsol经典实例012: 高斯波速的二次谐波产生激光系统是现代电子技术中的一个重要应用领域。

激光束的生成方法有很多种, 这些方法有个共同点: 波长由受激发射决定, 而受激发射取决于材料参数。

通常很难生成具有短波长的激光。

但是, 如果使用非线性材料, 就有可能产生频率为激光频率数倍的谐波。

通过使用二阶非线性材料可生成波长为基频光束波长一半的相干光。

本案例演示了如何设置非线性材料属性, 通过瞬态波仿真产生二次谐波。

模型中一束波长为1.06μm的激光聚焦于非线性晶体, 光束的腰部落于晶体内。

在激光的传播过程中, 大部分能量都集中在传播轴附近, 在求解麦克斯韦方程时可以近轴近似, 由此获得高斯波束。

一、物理场选择及预设研究Step01: 打开comsol软件, 单击“模型向导”选项创建模型, 在模型的“选择空间维度”界面选择“二维”, 在“选择物理场”界面分别选择“光学→波动光学→电磁波, 瞬态（ewt）”, 单击“添加”按钮。

对应变量设置完毕以后, 单击“研究”按钮, 在“选择研究”树中添加“预设研究”中的“瞬态”研究, 单击“完成”按钮进入软件主界面, 如图1所示。

图1 软件主界面二、全局定义1.参数Step02: 参数设置。

在模型开发器窗口的全局定义节点下, 单击“参数”子节点, 在“参数”设置窗口中, 定位到“参数”栏, 输入如图2所示的参数。

图2 设置全局参数2.解析定义Step03: 在“主屏幕”工具栏中单击“函数”选项, 在下拉菜单中选择“全局→解析”选项。

单击“解析1”子节点, 在“解析”设置窗口中, 定位到“函数名称”栏, 在文本输入框中输入“w”；定位到“定义”栏, 在“表达式”文本输入框中输入“w0*sqrt(1+(x/x)^2)”；定位到“单位”栏, 在“变元”文本输入框中输入“m”, 在“函数”在文本输入框中输入“m”,如图3所示。

Step04：在“主屏幕”工具栏中单击“函数”选项, 在下拉菜单中选择“全局→解析”选项。

COMSOL User Tips for Optical SimulationsHonghui ShenPhotonics Research Group(INTEC),Ghent University-IMEC,Sint-Pietersnieuwstraat41,9000Gent,Belgiumhonghuishen@July,20121MeshUsually the mesh size is set to beλ/6,where theλis the wavelength in the material. However for plasmonics the mesh size in the vicinity of the metal surface has to befiner, since thefield is confined to the metal surface by exciting the surface ually for silver in the visible wavelength range below5nm mesh size has to be applied to ensure the accuracy of simulations.One can do a convergency test for the mesh size around metal to avoid toofine or coarse meshes used.Since toofine mesh will definitely increase the calculation time and memory demand,on the other hand too coarse mesh will give non-accuratefield result around metal surface.For plasmonics,one has to avoid sharp corners of metal byfillet.Sharp corner will increase the mesh number and cause singularity.2Materials propertiesFor a simulation with wavelength sweep and dispersive material,the optical constants (or the permittivities)of materials corresponding to each wavelength has to be input into the model.There are several ways to do this:ing a analytical expression for the optical constant if ing a Drude-Lorentz dispersion model for metal tofit the measured optical constant1,for other materials like silicon the dispersive formula can be found in some literatures or online();1In COMSOL version4.3,this is implemented directly,see example Nanorods in v4.3.1ing the measured optical constant directly.For thefirst method,which is usually used in FDTD,the unknown parameters in Drude model have to be solved byfitting to the measured data.At some wavelength ranges the derivation of thefitting data by the Drude model is large with respect to the measured SOL provides a more easier way by using the measured data directly,i.e.the second method.For the second method,we have to build two separated txtfiles for both the real and the imaginary parts of the optical constant of e.g.silver.Be aware that in COMSOL a harmonic time dependence of e jωt is used,therefore for a loss material the optical constant is given by n-jk,where k>0.Figure1shows the txtfile for the real part of the optical constant of silver,n.Similarfile for k has to be built too.Figure1:Thefile example of optical constant of silver(dispersive material).Thefirst column is the wavelength in unit m,while the second column is the corresponding real part of refractive index of silver.For COMSOL v3.5,follow the steps below to input these twofiles(n,kfiles of optical constant)into model and set the corresponding domain:1.Clicking‘Option’→‘Functions...’,open the‘Functions’dialog box;2.Clicking the button‘New...’,open the‘New Function’dialog box(Figure2);3.Type e.g.n ag(the real part of optical constant of Ag,the function name can berandom)in the Function name editfield.Choose‘Interpolation’,and choose‘file’in the‘Use data from’(Figure2).4.Click‘Browse’,go to the directory where the nfile locates,import thefile;5.Following the above step,build another function for k ag(the imaginary part ofoptical constant of Ag);26.To this step two function n ag and k ag have been created.Now go to the sub-domain mode and open the’Subdomain settings’dialog box.Follow thefigure3 type n ag(lambda u)−j∗k ag(lambda u)in‘n’.Figure2:Input the optical constant into model as functionsFor COMSOL version above4,following the steps below2:1.Right click Materials node,create Material1(mat1)by choosing‘Material’.Rename‘Material1’as‘Ag’.2.Select the objects that correspond to material silver,choose Ag(mat1)andlocate and expand Material Properties in the Material Browser window.Select the‘Refractive index’under the‘Electromagnetic Models’and click the add sign below.Then a Refractive index(rfi)sub-node is added to the Ag (mat1)node.3.Right click Refractive index(rfi)and select‘Functions’→‘Interpolation’.AddInterpolation1(int1)to the node Refractive index(rfi).4.In the Interpolation Browser window,locate the parameters.Type‘n ag’inthe‘Function name’editfield,then choose‘file’in the‘Data source’.Click the ‘Browse...’button,go to the directory where the nfile locates and import the file,click‘open’to input thefile path into the‘Filename’editfield,and then click ‘Import’button.After this a table of real part of Ag data will be shown.2here only give one way to define material properties,there are other ways.3Figure3:Set the corresponding domain with the predefined functions.5.Repeat the last two steps for the imaginary part.6.After defining all of the materials used in current model in the Materials node,there is one more step to activate the material defined.In the Electromagnetic waves node,locate Material Type and select‘From material’,locate Electric Displacement Field and select‘Refractive index’in the‘Electric displacement field model’.3Light sourceThere are several different types of light source,e.g.point source,line source,plane source.Here ways to define plane source will be introduced.There are two usually used plane source,plane wave source and gaussian beam.For objects in a homogeneous material like a nanoparticle in air,plane wave can be introduced by using scattered formula,more details see examples‘Radar Cross Section’for COMSOL version above 4,and‘Dielectric scattering PML’for COMSOL v3.5.For gaussian beam,see examples‘Second Harmonic Generation of a Gaussian Beam’and‘nanorods’for COMSOL verson4.3,‘Second Harmonic Generation of a Gaussian Beam’for COMSOL v3.5.For models with object on a substrate,e.g.nanoparticles on glass,in these case,the scattered formula does not valid any more(since the object is not enclosed in a single4material).Several possible methods can be used:e two studies(see model example,‘plasmonic wire grating’in COMSOL v4.x),but the drawback of this method is that multiple ports have to be defined for each scattering order,it is feasible for cases with only few scattering orders existing.Ifa model has a lot of scattering orders,it becomes complicated to define all of theports for different orders.e the‘assembly’to introduce a soft source,see‘Assembly’in the‘Examples’fold(Details can be found in the slides included in‘Assembly’.This is also an example for how to set periodic boundary condition).But this seems only valid for COMSOL v3.5.Maybe this is implemented in the newest version4.3.3.Still use the scattered formula with adjustment applied to the initialfield(corre-sponding to the‘Incidentfield’in COMSOL v3.5,and the‘backgroundfield’in COMSOL v4.x).See‘Adjusted scattered’in the‘Examples’fold.4.(For other new methods,please add it here...)Now here details about method3in COMSOL v4.x(for v3.5it is similar)will be introduced.For simplicity,wefirst consider a2D model as shown infigure4(a), a periodic grating(material3,only one period is used in model therefore periodic boundary condition has applied to the lateral sides)on substrate(material2).In order to use the scattered formula,we have to givefield distribution in structure without grating,i.e.figure4(b)as the initialfiled of the model in(a).Thefield distribution in(b)can be easily calculated by Fresnel’s law together with Snell’s law.To simplify further we consider normal incidence.Therefore the totalfield in(b)can be expressed as:total=inc(y>y0)+refl(y>y0)+trans(y≤y0),(1)where‘y0’is the coordinate of the surface of substrate,‘inc’is the incident plane wave,‘refl’the reflected plane wave and the‘trans’the transmitted plane wave.Mathemat-ically Eq.2is realized with function‘sign’by:total=(inc+refl)·u1+trans·u2,(2)u1=(1+sign(y−y0))/2(3)u2=(1−sign(y+y0))/2(4)5(a)(b)Figure4:(a)A grating on a substrate.(b)structure without grating.These variables in the above3equations can be defined in the Variables under the Global Definitions node as shown infigure5.After define the initialfield expression in Variables,the initial(or background)field expression has to befilled into the Backgroundfiled as shown infigure6.More details see‘Adjusted scattered’in the ‘Examples’fold.4ScriptThe best way to learn script is saving an simple model as a script,i.e.the‘.m’file, which run on COMSOL with MATLAB.In a typical script,it usually includes the ge-ometry,mesh,physical setting(point(2D,3D),edge(3D),boundary(2D,3D)and sub-domain(2D,3D)setting),running,and postprocessing.In the following tips in terms of these aspects for COMSOL v3.5will be introduced(script for v4.x is quite different from v3.5and will be not introduced here,one is referred to the COMSOL documentation).4.1Create geometryGeometry example(the whole script is located:‘Scripts\Grating scattered formula’): %Geometrypml1=rect2(period,t pml);–create PML layer,rectangle6Figure5:Define initialfield in the Variables node.Figure6:Define the backgroundfield in the Backgroundfiled.7sub=rect2(period,t sub,‘base’,‘corner’,‘pos’,{‘0’,pos sub});–create PML refl=rect2(period,t refl,‘base’,‘corner’,‘pos’,{‘0’,pos refl});–create reflector p3ht=rect2(period,t p3ht,‘base’,‘corner’,‘pos’,{‘0’,pos p3ht});–create active layersup=rect2(period,t sup,‘base’,‘corner’,‘pos’,{‘0’,pos sup});–create air layer above active layerpml2=rect2(period,t pml,‘base’,‘corner’,‘pos’,{‘0’,pos pml});–create PMLg7=rect2(w gf/2,t gf,‘base’,‘corner’,‘pos’,{‘0’,pos gf y});–create grating on active layerfg=fillet(g7,‘radii’,ffillet,‘point’,[2,3]);–fillet grating sharp corners to avoid singularities.all=sub+refl+p3ht+sup+fg+pml1+pml2;–combine all the geometries created into one.%Geometry objects(following lines are usually needed in a COMSOL script to create geometry.Variable‘fem’includes all of the data in the model,e.g.geometry, physical setting,simulation results.)clear ss.objs={all};={‘all’};s.tags={‘all’};fem.draw=struct(‘s’,s);fem.geom=geomcsg(fem);figure,geomplot(fem);–use to visualize the geometry,to monitor if everything is correct.Tips:1.The‘base’,‘corner’in the rect2command can be neglected.The‘pos’can bespecified by numeric matrix([0,pos refl])instead of cell(e.g.{‘0’,pos refl});2.More command details can be found in the MATLAB command window by type‘help command name’;e command‘fillet’to round the sharp corners of metal components.mand‘geomplot’can be used to visualize the point,edge,boundary andsubdomain labels,which is useful when set corresponding physical conditions.8(Use the property‘edgelabel’of‘geomplot’to show edge label(edge in2D isboundary)).4.2Physical settingTo set conditions for subdomains,boundarys,edges,and points in a script,it is necessaryto label each of these items to distinguish them.The methods for setting conditions for these items are similar,so here we take the boundary and subdomain as examples(see example in‘Scripts\Grating scattered formula’).Before setting the boundary conditions,we have to know the boundary labels,andthis can be easily done by using‘geomplot’with the property‘edgelabel’:figure,geomplot(fem,‘edgelabel’,‘on’);–use to visualize the geometry togetherwith boundary labels.Or you can count it yourself there is a rule how the COMSOL number boundaries, counting from the low-left side to the up-right,but straight linefirst,then curves.After knowing the boundary labels,as following we can set the boundary conditions.clear bnd–clear the‘bnd’which is used to set the boundary condition.bnd.type={‘E0’,‘cont’,‘SC’};–specify in a cell how many different type of bound-ary conditions used in the model.‘E0’represents perfect boundary condition,‘cont’for continuity boundary condition,‘SC’for scattering boundary condition.(‘H0’for perfect boundary condition,but not used in the model.)For each type of the boundary condition in the cell has its own index,e.g.‘1’for‘E0’,‘2’for‘cont’...bnd.ind=[1,3,1,2,1,2,1,2,1,2,1,2,1,2,3,2,2,1,1,1,1,1,1,2,2];–specify the boundary conditions for each of the boundary in the geometry using the boundary condition index specify in the‘bnd.type’.As seen they are25boundaries,so we have to specify each of the boundaries in the sequence of boundary labels.appl.bnd=bnd;–boundary condition settingfinished.For subdomain condition,it is set by‘appl.equ’as following,see also in the example script.clear equequ.Sd={{‘Sdx guess rfweh’;‘lambdaS rfweh’},{‘Sdx guess rfweh’;‘Sdy guess rfweh’},9{‘Sdx guess rfweh’;‘Sdy guess rfweh’},{‘Sdx guess rfweh’;‘Sdy guess rfweh’}};–sub-domain setting,this is PML related,specify the width of PML in x-,y-axis directions.‘Sdx guess rfweh’,‘Sdy guess rfweh’and‘lambdaS rfweh’are some predefined value in COMSOL.For thefirst two means the PML widths in each direction use the geometric width.‘lambdaS rfweh’is used in the case when the geometric width is comparable with wavelength.This command line together with the next command line’equ.coordOn’are only used when there is a PML used in the model The number of types has to be the same number as the types of subdomain conditins specified in‘equ.n’equ.coordOn={{0;1},{0;0},{0;0},{0;0}};–specify which type of PML is activated in this model.this means thefirst‘n air’(stands for air)iin the‘equ.n’is set as PML. The‘{0;1}’stands for the PML absorb light propagating in the y-axis direction.The ‘0;0’with two‘0’deactivates the PML property in subdomain.equ.n={‘n air’,‘n air’,‘nag’,‘np3ht’};–specify subdomain material,‘n air’is defined in constants.‘n ag’and‘np3ht’are defined in global expressions;equ.Stype={‘coord’,‘none’,‘none’,‘none’};–PML type setting,the number of types has to be the same number as the types of subdomain conditins specified in‘equ.n’.‘coord’represents cartesian type of PML used,‘none’for no PML used.equ.matparams=‘n’;–subdomain material parameters,specify by refractive index ‘n’or permittivity.equ.ind=[1,2,3,4,3,2,1];–set the subdomain.appl.equ=equ;–subdomain settingfinished.4.3Geometric sweepFor some simulations geometric parameters have to been scanned to get an optimal design.For this can be accomplished by adding the Parametric Sweep into the Study node in COMSOL v4.x.But for v3.5,one can use COMSOL script(COMSOL with MATLAB)to sweep the geometric paramters with loops,see example for geometric sweep in‘Scripts\Grating scattered formula’.4.4PostprocessingCOMSOL GUI already integrates some useful postprocessing.However for some special user-defined parameters one can export the data to MATLAB,then process fundamental data calculated by COMSOL using script.Some important commands(in COMSOL10v3.5,command names are different in v4.x,corresponding command see documentation in v4.x)usually used in script are listed in following:Command(v3.5)FunctionGeomplot visualize the geometry of the model.Meshplot visualize the mesh.Postint integration in a boundary or domain.Postinterp extract the variable value e.g.E and H at specified points.Details about these commands type‘help command name’in the command window in COMSOL with Matlab.5Cluster runningMore information about high performance computing infrastructure at UGent see http: //www.ugent.be/hpc/en/infrastructure.Details on how to run COMSOL on cluster see user tips at userwiki http://hpc. ugent.be/userwiki/index.php/Main_Page or below.First we build the COMSOL model e.g.‘test22.mph’,then a job script is written to submit the job(the model will run on cluster)to the cluster.Notes on how to write a job script see userwiki.The job will wait in a queue before its running.Note that only job scripts can be submitted to the queuing system.These scripts can be written in a variety of languages e.g.bash,perl,python,...Following are the scripts used to submit jobs for running on single node(for model with low memory consumption,<16GB) and multinode with MPI(for model with large memory consumption,>16GB).COMSOL single node#!/bin/bash#PBS-l nodes=1:ppn=8–Use1node,8processes.#PBS-l walltime=00:30:00–Define the maximum walltime,30minutes,for the job to run.By default,this is1hour.cd$PBS O WORKDIRmodule load COMSOL/4.1.0.112–load COMSOLcomsol batch-np8-inputfile test22.mph–run the model‘test22.mph’COMSOL multinode job with MPI11There have been made several changes to the COMSOL start script.To start MPI jobs correctly,please use the following script example:#!/bin/bash#PBS-l nodes=3:ppn=8–Use3nodes and8processes on each node.#PBS-l walltime=00:30:00cd$PBS O WORKDIRmodule load ictce/4.0.6module load scriptsmodule load COMSOL/4.1.0.112##important.don’t change these variablesexport MYMPICMD=impirunexport MYMPIRUN HYBRID=1comsol-clustersimple batch-inputfile test22.mphAs for geometric sweep on cluster,MATLAB script should also work(but having no experience on this).To do a sweep with GUI,for v3.5,we have to build models separately for different geometric parameters,this can be simply done by saving‘fem’structure as mphfile with script command:flsave(‘modelname.mph’,fem)in a script to avoid building mphfile from GUI step by step.For version4.x,the geometric parameter sweep is embedded in the GUI directly,so one model is enough for different geometric parameters.But for a large model,it is better to build separate mph model for each geometric parameter.Because large model is already take quite long time to calculate,with geometric sweep more time needed.This means longer walltime needed,resulting in longer waiting time in the queue.In addition,longer time perhaps will increase the possibility of crash in the middle of running,resulting in losing data for these sweep parameters which have alreadyfinished.12。