光子晶體仿真大概是電磁學(xué)里“入門簡(jiǎn)單、精通極難”的典型。不少人第一次打開COMSOL覺(jué)得幾何能畫、邊界條件能設(shè)問(wèn)題不大真正上手做能帶才知道布里淵區(qū)掃描、布洛赫周期條件、特征頻率排序每一步都能卡掉一大批人。最近我把手頭這套基于COMSOL 5.6整理的光子晶體案例復(fù)現(xiàn)工程完整梳理了一遍——40多個(gè)mph文件覆蓋一維多層膜、二維晶格、三維木堆和反蛋白石結(jié)構(gòu)全部對(duì)應(yīng)那本常翻的光子晶體專著中的算例。這篇文章就想把里面的案例分布、參數(shù)邏輯和踩坑記錄都攤開講一講重點(diǎn)聊一聊用弱形式方程求解色散光子晶體能帶的思路希望對(duì)正在啃仿真、畫能帶、調(diào)網(wǎng)格的你有點(diǎn)實(shí)際幫助。1. 項(xiàng)目定位與整體價(jià)值拆解1.1 這40多個(gè)mph文件里的案例分布第一次拿到這套文件時(shí)我的第一反應(yīng)是“好東西終于不用從零開始搭幾何了”。逐個(gè)打開看了一遍之后才發(fā)現(xiàn)它其實(shí)有很清晰的梯度一維的六七個(gè)主要圍繞布拉格反射、一維禁帶、介質(zhì)膜堆透射譜二維的接近二十個(gè)正方晶格、三角晶格、Kagome晶格都有核心算例集中在色散關(guān)系、投影能帶、線缺陷波導(dǎo)、點(diǎn)缺陷微腔三維的十來(lái)個(gè)包括鉆石結(jié)構(gòu)、蛋白石/反蛋白石和木堆結(jié)構(gòu)。每個(gè)mph文件都自帶完整的幾何、物理場(chǎng)、網(wǎng)格和研究序列打開就能直接復(fù)現(xiàn)書上的結(jié)果。這個(gè)分布是有講究的。一維案例數(shù)量最少因?yàn)樗鼈儽举|(zhì)上是入門熱身物理過(guò)程最直觀計(jì)算量也最小主要用來(lái)建立“周期性結(jié)構(gòu)產(chǎn)生帶隙”的基本直覺(jué)。二維案例數(shù)量最多因?yàn)槎S是光子晶體研究的主戰(zhàn)場(chǎng)——既能講清楚能帶理論又不會(huì)像三維那樣動(dòng)不動(dòng)就吃掉幾十GB內(nèi)存。三維案例放在最后每個(gè)都代表一類經(jīng)典三維光子帶隙結(jié)構(gòu)計(jì)算代價(jià)高、收斂難但也最接近實(shí)際器件的真實(shí)情況。我在整理過(guò)程中發(fā)現(xiàn)文件名本身也透露了不少信息。比如一維的文件通常叫“Bragg_Mirror_Reflectance”或者“1D_PC_band”二維的會(huì)區(qū)分“TM”和“TE”偏振三維的則以結(jié)構(gòu)類型命名像“Woodpile_band”和“Inverse_opal”。如果你也是一邊看文獻(xiàn)一邊復(fù)現(xiàn)建議養(yǎng)成類似的文件命名習(xí)慣否則案例一多光靠記憶找模型會(huì)非常痛苦。1.2 為什么復(fù)現(xiàn)紙書案例比自建模型更有學(xué)習(xí)價(jià)值很多人喜歡從零開始建模型覺(jué)得這樣才叫“真本事”。我的看法相反在光子晶體這種理論、數(shù)值、物理圖像三者強(qiáng)耦合的領(lǐng)域復(fù)現(xiàn)經(jīng)典專著里的算例效率遠(yuǎn)高于自己悶頭造輪子。原因很簡(jiǎn)單紙書上的能帶圖和透射譜就是一份精確的“參考答案”。你自己搭的模型在參數(shù)設(shè)置、物理場(chǎng)選擇、邊界條件處理上有沒(méi)有問(wèn)題和書上的結(jié)果一比就暴露了。比如二維三角晶格光子晶體在空氣孔半徑r0.3a時(shí)TM禁帶出現(xiàn)在哪個(gè)頻率區(qū)間書上寫得清清楚楚。如果你算出來(lái)的帶隙位置偏移了10%以上那說(shuō)明要么填充率設(shè)錯(cuò)了要么波矢掃描路徑?jīng)]有按照高對(duì)稱點(diǎn)走要么網(wǎng)格粗到已經(jīng)影響精度。這種“對(duì)照糾錯(cuò)”的反饋速度是自己新建模型完全做不到的。mph文件還有一個(gè)紙面公式給不了的價(jià)值它把幾何、物理場(chǎng)、網(wǎng)格、求解器和后處理全部綁定在了一起?？凑撐哪憧赡苤恢馈安捎闷矫娌ㄕ归_法計(jì)算能帶”但具體到COMSOL里是用“電磁波頻域”接口做特征頻率還是用“弱形式PDE”寫方程網(wǎng)格單元是一階還是二階差分掃描步長(zhǎng)取多大這些細(xì)節(jié)才是決定結(jié)果能否復(fù)現(xiàn)的關(guān)鍵。而mph文件把這些“不可見的知識(shí)”完全暴露出來(lái)相當(dāng)于把一個(gè)熟練工程師的操作臺(tái)直接搬到你面前。1.3 從入門到進(jìn)階文件的階梯式用途這套案例庫(kù)對(duì)不同基礎(chǔ)的人價(jià)值點(diǎn)完全不一樣。剛接觸COMSOL的新手建議老老實(shí)實(shí)從一維案例開始先別急著改參數(shù)就把文件打開點(diǎn)“研究→計(jì)算”看看結(jié)果圖是怎么出來(lái)的。然后試著改一改介質(zhì)層的厚度觀察反射峰的位置怎么移動(dòng)。這個(gè)階段的目標(biāo)是熟悉軟件的基本操作流不要把精力花在理論上。有一定有限元基礎(chǔ)、但沒(méi)做過(guò)光子晶體的人重點(diǎn)應(yīng)該放在二維案例上。仔細(xì)研究模型樹里的“周期條件”節(jié)點(diǎn)和“輔助掃描”設(shè)置搞清楚波矢kx、ky和布里淵區(qū)路徑是怎么對(duì)應(yīng)的。能讀懂這些你就掌握了光子晶體能帶計(jì)算的主干流程。如果你已經(jīng)獨(dú)立算過(guò)幾個(gè)能帶想更進(jìn)一步那我建議你直接去啃三維案例和弱形式那部分。三維模型的網(wǎng)格剖分策略、特征頻率求解器的內(nèi)存管理、色散材料下非標(biāo)準(zhǔn)本征值問(wèn)題的處理這些是書上不會(huì)細(xì)講、但實(shí)際做科研幾乎必踩的技術(shù)點(diǎn)。把這層窗戶紙捅破你和“只會(huì)照看教程跑通模型”的人差距就拉開了。2. 一維、二維、三維建模的核心差異與方案選型2.1 一維物理場(chǎng)景最小門道卻不少一維光子晶體最常見的形式就是多層介質(zhì)膜高低折射率材料交替排列。這類結(jié)構(gòu)的仿真在COMSOL里實(shí)現(xiàn)起來(lái)非常直接幾何上畫一個(gè)矩形單元兩側(cè)設(shè)置周期性邊界條件上下或者左右設(shè)置端口邊界然后在頻域里掃一個(gè)寬頻段就能得到反射率和透射率譜。但“直接”不代表“簡(jiǎn)單”。我復(fù)現(xiàn)一維案例時(shí)第一個(gè)踩的坑就是端口邊界和周期性邊界條件的配合。計(jì)算透射譜時(shí)如果只在一側(cè)設(shè)置端口另一側(cè)直接設(shè)為地面得到的反射譜會(huì)和書上有明顯偏差因?yàn)樵谥芷诮Y(jié)構(gòu)里入射光會(huì)同時(shí)產(chǎn)生反射和透射能量守恒要求你在兩個(gè)端口上同時(shí)提取S參數(shù)。解決方法是使用“散射邊界條件”或顯式定義“端口”邊界把入射端口和出射端口都設(shè)置好后處理時(shí)用S11和S21分別表示反射率和透射率這樣復(fù)現(xiàn)出來(lái)的曲線才和書上一一對(duì)應(yīng)。一維的能帶計(jì)算又是另一種玩法。嚴(yán)格來(lái)說(shuō)一維周期結(jié)構(gòu)的能帶圖不需要完整的三維電磁場(chǎng)仿真可以把問(wèn)題退化成簡(jiǎn)單的特征頻率計(jì)算在單胞兩側(cè)加一對(duì)Floquet周期條件掃描波矢kx然后提取特征頻率。這種做法最見功夫的地方在于“歸一化頻率”的處理。紙書上習(xí)慣用ωa/2πc作為縱軸其中a是晶格常數(shù)ω是角頻率c是光速。你在COMSOL里算出來(lái)的特征頻率單位是Hz必須手動(dòng)除以c、乘以晶格常數(shù)a才能得到無(wú)量綱的歸一化頻率才能和文獻(xiàn)圖對(duì)上。這個(gè)換算關(guān)系我在早期復(fù)現(xiàn)時(shí)錯(cuò)了好幾次最后干脆寫了個(gè)參數(shù)“f_norm f * a / c_const”放在全局參數(shù)里一勞永逸。復(fù)用案例時(shí)還有個(gè)容易被忽略的點(diǎn)一維模型雖然物理維度低但COMSOL的“電磁波頻域”接口是三維求解器計(jì)算域仍然需要設(shè)置厚度。很多案例會(huì)把這個(gè)厚度設(shè)得很小比如0.01a本質(zhì)上是在用準(zhǔn)二維模型近似一維結(jié)構(gòu)。這種做法省內(nèi)存但如果你在研究?jī)A斜入射問(wèn)題就必須把厚度方向也設(shè)為周期條件否則就會(huì)出現(xiàn)人為的波導(dǎo)模式。2.2 二維是光子晶體仿真的主戰(zhàn)場(chǎng)二維案例是整個(gè)案例庫(kù)的核心也是我認(rèn)為最值得花時(shí)間吃透的部分。二維光子晶體的典型結(jié)構(gòu)是介質(zhì)柱正方晶格或者空氣孔三角晶格。在COMSOL里建模時(shí)只需要畫一個(gè)單胞比如正方形或者正六邊形把四邊設(shè)置成Floquet周期邊界條件然后在布里淵區(qū)路徑上掃描波矢就能得到完整的能帶圖。這里最關(guān)鍵的物理概念是偏振解耦。二維光子晶體中TM模式只有電場(chǎng)z分量Hz0TE模式只有磁場(chǎng)z分量Ez0兩者完全獨(dú)立。在做案例復(fù)現(xiàn)時(shí)你通常需要分別計(jì)算TM和TE能帶然后再把結(jié)果疊到同一張圖里看禁帶分布。COMSOL對(duì)此沒(méi)有自動(dòng)化的方法你需要手動(dòng)切換物理場(chǎng)設(shè)定——最常用的做法是定義兩個(gè)不同的“電磁波頻域”接口一個(gè)設(shè)置為“面外電場(chǎng)”TM另一個(gè)設(shè)置為“面外磁場(chǎng)”TE也可以分開跑兩個(gè)研究。部分案例的mph文件里直接把這兩種情況都做成了“組件”在“計(jì)算”前切換選擇即可復(fù)現(xiàn)時(shí)注意看清楚當(dāng)前是哪個(gè)偏振在起作用。二維能帶計(jì)算中布里淵區(qū)掃描路徑的設(shè)定對(duì)初學(xué)者來(lái)說(shuō)是最難理解的一環(huán)。以正方晶格為例第一布里淵區(qū)的高對(duì)稱點(diǎn)為Γ(0,0)、X(0.5,0)、M(0.5,0.5)歸一化坐標(biāo)能帶圖通常沿?！鶻→M→Γ這條折線掃。操作上你需要引入一個(gè)路徑參數(shù)把波矢分量kx、ky寫成路徑參數(shù)的函數(shù)再用“輔助掃描”去掃描它。我在復(fù)現(xiàn)時(shí)踩過(guò)的一個(gè)坑是有的案例在定義Floquet周期條件時(shí)會(huì)把kx和ky寫反導(dǎo)致能帶扭曲。排查方法很簡(jiǎn)單把某個(gè)高對(duì)稱點(diǎn)的頻率和書上對(duì)比如果Γ點(diǎn)處兩個(gè)不同方向的模式頻率不相等那說(shuō)明k矢量的分量映射出問(wèn)題了。二維案例里還有一大批是做缺陷態(tài)的。“線缺陷波導(dǎo)”和“點(diǎn)缺陷微腔”都需要把單胞擴(kuò)展成超胞比如5×5的介質(zhì)柱陣列中間去掉一根柱子或者一整行柱子然后再加一圈周期邊界條件。超胞模型的網(wǎng)格數(shù)量會(huì)迅速上升計(jì)算時(shí)間和內(nèi)存也相應(yīng)上漲。我復(fù)現(xiàn)線缺陷波導(dǎo)案例時(shí)最開始直接用了默認(rèn)網(wǎng)格結(jié)果算出來(lái)的導(dǎo)帶模式頻率偏高了接近8%后來(lái)把缺陷附近的網(wǎng)格手動(dòng)加密到λ/12左右才得到和書吻合的色散曲線。做這類案例時(shí)“局部網(wǎng)格加密”是非常關(guān)鍵的一步而且要學(xué)會(huì)利用COMSOL的“細(xì)化”工具不必全局加密只在缺陷區(qū)域細(xì)劃就能兼顧精度和速度。2.3 三維模型的一切難點(diǎn)都在“代價(jià)”上三維光子晶體案例比如金剛石結(jié)構(gòu)、木堆結(jié)構(gòu)和反蛋白石結(jié)構(gòu)是整套mph文件里最“勸退”的部分。不是說(shuō)物理上多復(fù)雜而是計(jì)算代價(jià)呈數(shù)量級(jí)上漲。一個(gè)典型的木堆結(jié)構(gòu)單胞如果網(wǎng)格用二階四面體單元數(shù)輕松突破百萬(wàn)量級(jí)特征頻率研究需要的內(nèi)存常常超過(guò)32GB計(jì)算時(shí)長(zhǎng)以小時(shí)計(jì)。所以三維案例的復(fù)現(xiàn)策略和二維完全不一樣。在二維里你可以直接上高精度網(wǎng)格但在三維里必須先“舍得”——先用極粗的網(wǎng)格跑一遍只求能帶趨勢(shì)對(duì)確定參數(shù)沒(méi)問(wèn)題之后再逐步加密網(wǎng)格觀察關(guān)鍵禁帶邊緣頻率是否還有明顯漂移。COMSOL里的“參數(shù)化掃描”和“網(wǎng)格序列”可以聯(lián)動(dòng)你可以把網(wǎng)格最大單元尺寸設(shè)成一個(gè)掃描參數(shù)先跑0.3a、再跑0.2a、0.15a看結(jié)果收斂情況。這個(gè)收斂性檢查是我在復(fù)現(xiàn)三維能帶時(shí)習(xí)慣性做的一步能替你省下大量不必要的全精度計(jì)算。另外一個(gè)省內(nèi)存的技巧是充分利用結(jié)構(gòu)對(duì)稱性。很多三維光子晶體結(jié)構(gòu)具有平移、旋轉(zhuǎn)和鏡面對(duì)稱性你可以把計(jì)算域從整個(gè)單胞縮減為1/2或1/4然后在剖切面上設(shè)置對(duì)稱條件或周期性條件。比如反蛋白石結(jié)構(gòu)沿對(duì)角線對(duì)稱用半胞計(jì)算理論上內(nèi)存占用幾乎減半。需要注意不是所有結(jié)構(gòu)都能隨便切對(duì)三維布里淵區(qū)路徑上的某個(gè)波矢點(diǎn)對(duì)稱條件的具體形式可能不同。這點(diǎn)要結(jié)合群論知識(shí)仔細(xì)核對(duì)否則你算出來(lái)的能帶會(huì)平白多出一些“禁帶”又或者漏掉某些本該存在的模式。三維能帶的布里淵區(qū)路徑也更長(zhǎng)例如面心立方結(jié)構(gòu)的路徑是Γ-X-U-L-Γ-W-K波矢掃描點(diǎn)數(shù)通常需要50到100個(gè)才夠平滑。每個(gè)波矢點(diǎn)都要做一次特征頻率求解這就意味著一次完整的能帶計(jì)算相當(dāng)于跑50到100次高負(fù)載仿真。我的建議是先只跑低對(duì)稱點(diǎn)線路的粗網(wǎng)格結(jié)果確認(rèn)沒(méi)有明顯物理異常后再一次性提交完整掃描并且把COMSOL的“作業(yè)序列”和并行計(jì)算配好否則中途一斷電或者一次參數(shù)填錯(cuò)就要從頭再來(lái)心態(tài)很容易崩。3. 實(shí)操過(guò)程從打開mph文件到跑出能帶圖3.1 打開模型后先做這四件事拿到一個(gè)mph文件不要急著點(diǎn)“計(jì)算”。我每次復(fù)現(xiàn)新案例都會(huì)先按固定順序檢查四樣?xùn)|西效率高出很多。第一打開“全局定義參數(shù)”節(jié)點(diǎn)看一眼a晶格常數(shù)、r半徑、填充率、epsilon_r介電常數(shù)這些核心參數(shù)是否存在單位是什么。很多案例復(fù)現(xiàn)失敗根源只是參數(shù)單位弄混了比如把半徑寫成微米級(jí)的數(shù)值卻沒(méi)有同步設(shè)置幾何單位。第二展開“研究1”節(jié)點(diǎn)查看是否已經(jīng)有“輔助掃描”或“參數(shù)化掃描”步驟。如果有點(diǎn)開看掃描的變量名和取值列表。比如二維正方晶格案例通常會(huì)有一個(gè)名為“path”的路徑參數(shù)或者直接是變量“kx”。搞清楚掃描設(shè)置是理解整個(gè)模型如何生成能帶圖的關(guān)鍵。第三在“組件→定義”里找到“周期條件”或“周期性邊界條件”確認(rèn)邊界條件類型是“Floquet周期”并且檢查k矢量分量表達(dá)式。這里最容易出問(wèn)題因?yàn)槊字茊挝幌虏ㄊ竼挝皇莚ad/m而書上習(xí)慣用與晶格常數(shù)相關(guān)的歸一化單位你可能需要自己做換算。第四看網(wǎng)格的大概尺寸。在“網(wǎng)格→大小”里查看“最大單元尺寸”如果它是全局統(tǒng)一的就要特別注意缺陷案例是否局部加密。一維和二維案例對(duì)網(wǎng)格不太敏感三維案例對(duì)網(wǎng)格非常敏感先估算一下這個(gè)模型大概要多久才能算完避免盲目點(diǎn)下“計(jì)算”之后傻等。這四步檢查做完再點(diǎn)“研究→計(jì)算”通常就能順利復(fù)現(xiàn)。如果哪一次結(jié)果不對(duì)回查這四步往往能快速定位問(wèn)題。3.2 布洛赫周期條件與布里淵區(qū)掃描怎么配合布洛赫周期條件對(duì)剛接觸的人而言是光子晶體仿真里最“玄學(xué)”的一環(huán)。一句話解釋它告訴求解器單胞邊界上的場(chǎng)從一個(gè)邊界到另一個(gè)邊界只相差一個(gè)相位因子e^{ik·a}其中a是晶格矢量k就是你正在掃描的波矢。這個(gè)相位因子的存在把無(wú)限大的周期結(jié)構(gòu)問(wèn)題歸結(jié)到單個(gè)原胞上求解。在COMSOL的“電磁波頻域”模塊里添加這個(gè)條件不需要手動(dòng)寫任何復(fù)數(shù)表達(dá)式只需要在“周期性條件”節(jié)點(diǎn)類型中選擇“Floquet周期”然后設(shè)置k矢量的分量。以二維正方晶格為例如果晶格常數(shù)是a高對(duì)稱路徑上的點(diǎn)用歸一化坐標(biāo)(kx0, ky0)表示那么實(shí)際波矢分量為kx kx0 * 2π / aky ky0 * 2π / a比如X點(diǎn)歸一化坐標(biāo)是(0.5, 0)M點(diǎn)是(0.5, 0.5)Γ點(diǎn)是(0,0)。在參數(shù)化掃描中你通常定義一個(gè)路徑參數(shù)s讓s從0變到1然后讓(kx0, ky0)沿著?！鶻→M→Γ折線變化。具體實(shí)現(xiàn)時(shí)可以用條件表達(dá)式分段定義kx0 if(s0.4, s/0.40.5, if(s0.8, 0.5-(s-0.4)/0.40.5, (s-0.8)/0.20.5))ky0 if(s0.4, 0, if(s0.8, (s-0.4)/0.40.5, 0.5-(s-0.8)/0.2*0.5))雖然寫起來(lái)有點(diǎn)繞但這是能帶圖橫軸坐標(biāo)唯一的正確控制方式。還有一種稍省事的辦法直接定義“顯式”參數(shù)列表只用十幾個(gè)關(guān)鍵高對(duì)稱點(diǎn)每個(gè)點(diǎn)之間斜率不同圖形也能看但曲線不平滑禁帶邊界不夠清晰。實(shí)際復(fù)現(xiàn)時(shí)我傾向于至少掃描40~60個(gè)點(diǎn)這樣能帶上的拐點(diǎn)、簡(jiǎn)并點(diǎn)都能看得很準(zhǔn)。3.3 特征頻率研究排序、導(dǎo)出、畫圖特征頻率研究是能帶計(jì)算的核心。設(shè)置上建議在“特征頻率”研究的“設(shè)置”窗口中指定“搜索特征頻率的基準(zhǔn)值”和“目標(biāo)特征頻率數(shù)”?；鶞?zhǔn)值的選取很講究如果你只關(guān)心禁帶附近幾條帶可以把基準(zhǔn)值設(shè)為某個(gè)中心頻率讓求解器只找鄰域里的模式。如果你想要一張完整的能帶圖就要把目標(biāo)特征頻率數(shù)設(shè)多比如每點(diǎn)算20個(gè)模式保證高頻段的能帶也完整。實(shí)際運(yùn)行中最大的痛點(diǎn)不是算不出特征頻率而是算出來(lái)的模式不知道怎么對(duì)應(yīng)到一條連續(xù)的能帶上。特征頻率求解器在每個(gè)波矢點(diǎn)都獨(dú)立求解出來(lái)的頻率是從低到高排的但物理上應(yīng)該連續(xù)的能帶在數(shù)值上會(huì)因?yàn)槟Ｊ浇徊妗?duì)稱性簡(jiǎn)并而“斷掉”或者“跳帶”。比如在第5個(gè)波矢點(diǎn)第6和第7條帶可能是簡(jiǎn)并的到下一點(diǎn)它們可能徹底分開。如果你只是機(jī)械地把每個(gè)點(diǎn)的前6個(gè)頻率連成一條線畫出來(lái)的曲線會(huì)出現(xiàn)突跳怎么看都不對(duì)。解決這個(gè)問(wèn)題沒(méi)有萬(wàn)能公式但我有一個(gè)比較實(shí)用的組合拳先在COMSOL里把每個(gè)波矢點(diǎn)的所有特征頻率導(dǎo)出然后在繪圖階段按“模式序號(hào)”分組顯示觀察哪些模式具有相同的場(chǎng)分布特征比如位移場(chǎng)奇偶性、節(jié)線數(shù)人工把它們排到同一條能帶上。某些mph文件里會(huì)直接使用COMSOL的結(jié)果“全局計(jì)算”節(jié)點(diǎn)在設(shè)置里選中“特征頻率”并啟用“按模式追蹤”軟件會(huì)自動(dòng)跟蹤相鄰波矢點(diǎn)的連續(xù)模式。如果自動(dòng)追蹤失效還有一個(gè)土辦法——把波矢掃描步長(zhǎng)加密一倍看能否恢復(fù)連續(xù)通常跳變是因?yàn)閮刹街g的模式交換太劇烈加密后情況會(huì)明顯好轉(zhuǎn)。畫出能帶圖后橫軸是路徑參數(shù)s或者實(shí)際波矢距離縱軸是歸一化頻率。我習(xí)慣在COMSOL里用“一維繪圖組→點(diǎn)圖”把每個(gè)波矢點(diǎn)的多個(gè)特征頻率畫成離散點(diǎn)再用第三方軟件擬合連線。這樣能看出哪些點(diǎn)是簡(jiǎn)并的、哪些點(diǎn)是真正的禁帶邊緣信息量比直接連線大得多。3.4 弱形式方程求解色散光子晶體能帶進(jìn)階這部分是整套復(fù)現(xiàn)工程里我認(rèn)為最有“進(jìn)階感”的內(nèi)容也是最近很多人開始問(wèn)的東西——基于COMSOL弱形式方程求解色散光子晶體能帶。為什么要用弱形式因?yàn)闃?biāo)準(zhǔn)的“電磁波頻域”接口在求解特征頻率時(shí)默認(rèn)把介電常數(shù)看成一個(gè)常數(shù)本征值方程是線性的。但實(shí)際材料往往是色散的比如金屬在可見光和近紅外波段遵循Drude模型介電常數(shù)是頻率的函數(shù)ε(ω)。這種情況下本征方程變成了非線性本征值問(wèn)題直接用ewfd的“特征頻率”研究算出來(lái)的頻率根本不含材料色散或者說(shuō)只能在某個(gè)固定頻率處對(duì)ε取值做不到能帶上每一個(gè)k點(diǎn)都對(duì)應(yīng)正確頻率下的ε。弱形式可以繞開這個(gè)限制。思路是把布洛赫解直接“手動(dòng)”寫進(jìn)方程。電場(chǎng)寫成E(r)U(r)e^{ik·r}其中U是周期函數(shù)k是波矢。把這個(gè)形式代進(jìn)頻域亥姆霍茲方程后原來(lái)的旋度算子要替換成帶修正項(xiàng)的算子也就是把每個(gè)空間導(dǎo)數(shù)?變換為?ik。在弱形式PDE中你可以逐項(xiàng)寫出修正后的方程并把介電常數(shù)寫成和本征值λ即頻率ω相關(guān)的表達(dá)式讓求解器在每次迭代中都“重新考慮”材料參數(shù)。COMSOL里的具體操作流程是這樣的第一步添加“數(shù)學(xué)”→“PDE接口”→“弱形式PDE”。定義因變量為電場(chǎng)分量ex、ey、ez注意它們都是復(fù)數(shù)。第二步在“全局定義→變量”里定義波矢分量kx、ky以及材料參數(shù)例如Drude模型ε1-ωp2/(ω2iγω)。第三步在弱表達(dá)式中寫入積分方程核心是電場(chǎng)雙旋度項(xiàng)與介電項(xiàng)。第四步把研究類型設(shè)為“特征值”并把本征值變量λ映射成-ω2/c2或直接定義ω2-λ*c2這樣每解出一個(gè)λ就得到一個(gè)對(duì)應(yīng)的頻率。第五步設(shè)定波矢掃描路徑和前面ewfd的方式完全一致。這個(gè)流程我第一次走得并不順最大的問(wèn)題出在“特征值變換”上。COMSOL弱形式默認(rèn)的特征值量綱和你定義的物理量不一定對(duì)得上極容易算出量綱離譜的頻率。我的建議是在“研究→步驟→特征值→設(shè)置”里利用“特征值變換”把本征值λ的實(shí)部解釋為與頻率平方相關(guān)的數(shù)例如omega sqrt(-real(lambda)) * c / a然后把這些派生變量加到“全局計(jì)算”里輸出歸一化頻率直接畫圖。如果你發(fā)現(xiàn)算出來(lái)的能帶在Γ點(diǎn)附近頻率不為零或者為負(fù)多半是符號(hào)設(shè)置問(wèn)題把方程的符號(hào)整體換一下很快就能修正。用弱形式的好處不只是能處理色散材料還在于它讓你對(duì)“求解過(guò)程”本身有了完全的掌控。比如你想研究增益介質(zhì)、想加入磁光效應(yīng)、想改方程形式加入高階項(xiàng)標(biāo)準(zhǔn)物理場(chǎng)接口可能根本不給你這個(gè)機(jī)會(huì)但在弱形式下你只需要改動(dòng)一個(gè)積分表達(dá)式即可。壞處則是你必須自己負(fù)責(zé)正確的量綱和邊界條件調(diào)試門檻高。作為對(duì)比我的經(jīng)驗(yàn)是只要能直接用ewfd解決的非色散問(wèn)題就不要強(qiáng)行用弱形式一旦材料色散成為繞不開的坎弱形式就是值得投入的方向。這套案例庫(kù)里的弱形式文件跑通一遍并和書上的色散能帶對(duì)照過(guò)會(huì)對(duì)本構(gòu)關(guān)系與數(shù)值求解的關(guān)系有一個(gè)比讀十篇論文都深的體會(huì)。4. 常見問(wèn)題與排查技巧實(shí)錄4.1 能帶斷帶/跳變從哪里找原因能帶圖最輕的毛病是“斷帶”——相鄰兩個(gè)波矢點(diǎn)的頻率曲線不在同一個(gè)位置接上。絕大多數(shù)情況的根源是模式排序不一致。特征頻率求解器在每個(gè)波矢點(diǎn)都獨(dú)立求解并不會(huì)自動(dòng)識(shí)別“第n條帶”如果兩個(gè)波矢點(diǎn)之間模式發(fā)生交叉排序邏輯就會(huì)錯(cuò)位曲線看起來(lái)就像從一根跳到了另一根。面對(duì)斷帶問(wèn)題我一般按三層順序排查。第一層檢查掃描步長(zhǎng)是否太粗。路徑點(diǎn)太疏時(shí)高頻模式交叉頻繁跳變幾乎無(wú)法避免把掃描點(diǎn)翻倍往往能平滑很多。第二層檢查網(wǎng)格尺寸是否合理。粗網(wǎng)格導(dǎo)致低頻模式還算準(zhǔn)確但高頻模式嚴(yán)重畸變不同模式的誤差大小不一樣畫出來(lái)就差得更遠(yuǎn)。第三層如果你在結(jié)果里啟用了“按模式追蹤”檢查它的追蹤準(zhǔn)則是否和你的物理場(chǎng)景匹配。對(duì)于強(qiáng)耦合模式自動(dòng)追蹤會(huì)失敗這時(shí)只能手工分組處理。4.2 負(fù)頻率、虛數(shù)頻率和“幽靈?！彼闾卣黝l率時(shí)結(jié)果里偶爾會(huì)出現(xiàn)頻率為負(fù)或者虛部的數(shù)值第一次見到的確會(huì)慌但這類情況大多不是致命錯(cuò)誤。負(fù)頻率的出現(xiàn)通常只是特征值符號(hào)約定問(wèn)題。COMSOL在解本征方程時(shí)會(huì)同時(shí)解出正負(fù)對(duì)稱的特征值取絕對(duì)值或者確認(rèn)符號(hào)映射關(guān)系就行。真正需要警惕的是虛部明顯的頻率如果在損耗材料比如有損介質(zhì)、金屬中計(jì)算虛頻率反映的是該模式的實(shí)際損耗或增益這是物理信息但如果材料是無(wú)損的虛頻率卻仍然存在說(shuō)明網(wǎng)格精度不足產(chǎn)生了“幽靈?！?。鑒別“幽靈模”最有效的方法是網(wǎng)格收斂性檢查。把最大網(wǎng)格尺寸從λ/8加密到λ/16如果某個(gè)可疑頻率的實(shí)部變化超過(guò)1%而虛部沒(méi)有消失那基本可以認(rèn)為是偽解。還有一種高頻“幽靈?！眮?lái)自棱角奇異性比如介質(zhì)柱邊緣的場(chǎng)強(qiáng)趨于無(wú)限大數(shù)值上會(huì)擠出假的局域模式。處理方式是在幾何設(shè)計(jì)上加上微小的圓角或者局部加密奇異點(diǎn)周圍網(wǎng)格不要用非常銳利的棱邊去碰這類問(wèn)題。4.3 三維計(jì)算內(nèi)存不夠怎么硬扛三維案例復(fù)現(xiàn)時(shí)最難熬的是內(nèi)存。“Out of memory”這個(gè)紅色報(bào)錯(cuò)我在木堆結(jié)構(gòu)上見過(guò)很多次。COMSOL默認(rèn)的求解器通常是PARDISO或MUMPS直接分解稀疏矩陣雖然穩(wěn)但內(nèi)存開銷巨大。應(yīng)對(duì)策略按優(yōu)先級(jí)排序第一換迭代求解器。在“研究→求解器配置”中把“特征值求解器”改為“迭代”或“FEAST”。特征值問(wèn)題用FEAST算法往往比直接分解更快且內(nèi)存占用相對(duì)平緩。第二縮小目標(biāo)模式數(shù)。三維能帶不需要一次算20條帶剛開始先算最低的6~8條確認(rèn)能帶路徑正確后再逐步增加模式數(shù)。第三利用對(duì)稱性縮減模型尺寸。前面說(shuō)過(guò)的半胞、1/4胞策略在三維案例里能顯著降低網(wǎng)格量。第四如果你有計(jì)算集群或者幾個(gè)核心比較多的機(jī)器別忘了把COMSOL的“并行”打開至少跑起來(lái)不會(huì)卡死在單核上。另外有一個(gè)常被忽略的點(diǎn)網(wǎng)格的階次。三維四面體默認(rèn)常用二階拉格朗日單元精度好但代價(jià)高。如果內(nèi)存實(shí)在頂不住可以嘗試對(duì)遠(yuǎn)離結(jié)構(gòu)細(xì)節(jié)的區(qū)域使用一階單元或者直接在“網(wǎng)格”里使用“邊界層”方式對(duì)關(guān)鍵界面精細(xì)剖分其余部分粗剖。這樣算出來(lái)的能帶邊緣頻率和全二階模型差距通常小于5%但內(nèi)存需求可能下降一半。4.4 mph文件版本兼容與打開異常COMSOL的mph文件雖然自帶完整模型信息但不同主版本之間的兼容性并不完美。5.6版本的mph在6.x里一般能直接打開只是可能提示“求解器默認(rèn)設(shè)置將升級(jí)”反過(guò)來(lái)6.x保存的文件低版本打不開這點(diǎn)要特別注意。如果你手頭只有5.6別人發(fā)給你一個(gè)6.x的mph最省事的方法是讓對(duì)方在“另存為”時(shí)選擇“保存副本→版本5.6”。我在整理這套案例文件時(shí)還遇到過(guò)一種奇怪情況文件雙擊后打開模型樹是空的或者提示“需要一個(gè)或多個(gè)組件缺失”。這通常不是文件損壞而是模型引用了自定義材料庫(kù)或外部CAD幾何文件路徑失效了。解決辦法是重建缺失組件——把材料參數(shù)手動(dòng)敲進(jìn)去或者重新導(dǎo)入幾何文件。如果你準(zhǔn)備長(zhǎng)期保存這套復(fù)現(xiàn)工程強(qiáng)烈建議把每個(gè)案例用到的幾何文件和自定義材料定義都放在同一個(gè)文件夾里mph與外部依賴分離存儲(chǔ)這樣無(wú)論拷到哪臺(tái)機(jī)器上都能直接復(fù)現(xiàn)。最后補(bǔ)一個(gè)小提醒mph文件在計(jì)算過(guò)程中如果被強(qiáng)制終止比如斷電、任務(wù)管理器強(qiáng)殺很可能會(huì)因?yàn)閷懭氩煌暾鵁o(wú)法再次打開。養(yǎng)成“原文件只讀、復(fù)制工作副本再跑”的習(xí)慣可以避免很多心痛時(shí)刻。我在這些案例身上花的時(shí)間不短最深的體會(huì)是復(fù)現(xiàn)一套能帶圖真正值錢的不是那張彩色圖片本身而是你被迫去理解每一處參數(shù)、每一個(gè)邊界條件、每一次求解器選擇的過(guò)程。等你跑通了第三十個(gè)案例再拿到一個(gè)新結(jié)構(gòu)時(shí)你會(huì)自然地知道第一步應(yīng)該檢查什么、可能會(huì)死在哪個(gè)環(huán)節(jié)這種直覺(jué)是看多少篇教程都換不來(lái)的。如果后面你有機(jī)會(huì)折騰這套文件建議邊跑邊記筆記把每個(gè)案例的參數(shù)為什么要這樣設(shè)置、你踩了什么坑都寫下來(lái)這套整理過(guò)的知識(shí)比文件本身更值錢。