光子晶體仿真看起來門檻高實際上絕大多數(shù)時間都花在修正模型的細節(jié)上。我最早跟著論文復現(xiàn)二維空氣孔光子晶體整整一周都在跟能帶圖里的鋸齒較勁最后才發(fā)現(xiàn)只是材料介電常數(shù)虛部沒有清零。為了徹底擺脫反復試錯的局面我決定用COMSOL 5.6把《光子晶體》教材中的典型案例完整復現(xiàn)一遍。目前這套復現(xiàn)項目積累了40多個可直接運行的mph文件涵蓋一維、二維、三維三種維度的光子晶體結構包括透射譜、反射譜、能帶圖、本征模場分布等完整輸出。這篇文章既是對這套案例庫的說明也是把復現(xiàn)過程中踩過的坑、總結的方法和調(diào)參經(jīng)驗一次性講清楚。無論你是剛接觸光子晶體仿真的研究生還是已經(jīng)在算能帶但經(jīng)常對不上文獻結果的工程師這套從一維到三維的完整鏈路都值得參考。1. 為什么我要在COMSOL 5.6里復現(xiàn)光子晶體案例1.1 從“看書懂”到“動手會”的鴻溝光子晶體相關書籍通常會把能帶理論講得很細但真正動手建模的時候你會發(fā)現(xiàn)書上的信息根本不夠用。比如書里可能只寫了晶格常數(shù)a600 nm、空氣孔半徑r180 nm卻沒有寫清楚Floquet邊界條件的波矢量該怎么設置也沒有說明用TM模式還是TE模式計算。這些信息缺失導致一百個人能跑出一百種結果。所以我決定換一個思路不再零散地在網(wǎng)上找示例而是選一本案例最完整、參數(shù)最清晰的光子晶體專著把其中能復現(xiàn)的算例逐個做出來。所謂復現(xiàn)不是簡單畫出幾何而是讓計算得到的禁帶位置、能帶寬度、透射率曲線和書中的結果一致。有了這樣一條校準線后續(xù)做新結構的時候我才有底氣去改參數(shù)、換材料知道哪些結果是合理的哪些是模型出錯了。1.2 案例庫的組織方式與文件規(guī)范40多個mph文件如果隨意堆在一起半年后連自己都找不到。我在整理案例庫時采用“維度—結構類型—求解目標”的三級目錄結構根目錄下分1D、2D、3D三大類每一類再按照具體結構細分。比如2D目錄下就有正方晶格空氣孔、三角晶格空氣孔、六角晶格介質(zhì)柱、線缺陷波導、微腔諧振器等子目錄。文件命名也有一套固定的規(guī)則。我會把結構類型、關鍵幾何參數(shù)、折射率組合、計算模式四個要素寫進文件名。例如2D_Triangular_r0.3_n2.34_TM_band.mph看到名字就知道這是三角晶格空氣孔結構r/a0.3背景折射率2.34算的是TM模式能帶。另外每個子目錄里放一個README.md記錄案例來源、章節(jié)號、預期結果和模型注意事項。這樣即使隔幾個月再打開也能快速恢復上下文。1.3 選擇COMSOL 5.6的具體原因標題里的“Comsol56”指的就是COMSOL 5.6版本。我選擇這個版本主要是看中它的波光學模塊穩(wěn)定性。5.6對Floquet周期邊界、端口邊界和散射邊界的底層求解器做了不少優(yōu)化特征值計算的收斂性比早期版本強很多。特別是對高介電常數(shù)對比度的結構早期版本經(jīng)常出現(xiàn)“找不到模式”的提示5.6的容錯明顯變好。另一個重要原因是5.6的“數(shù)學模塊”里弱形式PDE接口增強了非線性求解能力。后面我要專門講到的“基于COMSOL弱形式方程求解色散光子晶體能帶”正是依賴這個接口。早期版本在這個功能上偏弱很多頻散材料模型需要額外寫代碼才能求解。5.6把弱形式的穩(wěn)定性提升了一個臺階才讓我能在原生界面里完成色散能帶計算。2. 光子晶體仿真的三條核心方法論2.1 布里淵區(qū)、倒格子與k空間路徑不管是幾維結構光子晶體仿真都繞不開三個核心概念正格子、倒格子、布里淵區(qū)。正格子是你在COMSOL里畫的周期性晶胞倒格子是與之對應的動量空間周期單元布里淵區(qū)則是倒格子的原胞。能帶圖描繪的就是本征頻率在這個布里淵區(qū)內(nèi)沿特定路徑的變化。初學者最容易混淆的是幾何尺寸和k點路徑之間的關系。幾何尺寸可以用實際晶格常數(shù)建模比如三角晶格a600 nm但能帶計算里的k點必須沿布里淵區(qū)邊界走。三角晶格的高對稱路徑是Γ(0,0) → M(0.5,0) → K(0.333,0.333) → Γ(0,0)這些坐標是無量綱的是相對于倒格子基矢的。如果直接把路徑坐標當作普通變量輸入COMSOL得到的能帶會整體變形。我的處理方法是把倒格子基矢的換算關系直接寫進模型的全局參數(shù)里。以三角晶格為例倒格子基矢長度為4π/(a√3)Floquet邊界需要的波矢量分量定義為kx 4*pi/(a*sqrt(3)) * s1 ky -4*pi/(3*a) * s1這樣掃描參數(shù)s1遍歷0到1的區(qū)間k點就自動沿高對稱路徑移動。我現(xiàn)在做二維案例時已經(jīng)把這套表達式做成了公共參數(shù)組復制到任意二維模型都能直接用只需根據(jù)晶格類型修改系數(shù)。2.2 Floquet周期邊界條件的參數(shù)化Floquet邊界條件是光子晶體能帶計算的基石。它的作用是讓晶胞兩側的電磁場滿足一個相位關系相當于把無限周期結構的邊界效應濃縮到一個單元里。在COMSOL中設置周期邊界時類型必須選“Floquet周期性”不能選普通的“周期性”否則邊界兩側的電場無法傳播相移。邊界條件界面里需要指定兩個方向的波矢量分量kFloq1和kFloq2。這兩個分量的單位是rad/m不是倒格子坐標。很多人在這一步出錯是因為直接把k點坐標輸入進去導致相位積累錯誤。正確做法是換算正方晶格中kFloq1 2π·kx/akFloq2 2π·ky/a三角晶格則需要考慮基矢夾角手動算出兩個方向的投影系數(shù)。在參數(shù)化掃描過程中我會把kFloq1和kFloq2定義成全局參數(shù)的表達式然后用“輔助掃描”功能讓k點連續(xù)遍歷高對稱路徑。這里有個經(jīng)驗掃描點數(shù)量不是越多越好。我通常設置61個掃描點既能保證能帶曲線平滑又不至于讓求解時間翻倍。對于多維參數(shù)掃描COMSOL的“參數(shù)掃描”會為每一組參數(shù)完整求解一次掃描點過多時建議拆成兩段執(zhí)行方便中途查看結果。2.3 特征值求解器的目標設置與模式篩選特征值求解器是能帶計算的核心引擎。COMSOL默認會計算“所需模式數(shù)”個最低頻率的特征模但光子晶體往往需要特定頻率區(qū)間的模式而不是最低的那幾個。我的習慣是把“特征值搜索范圍”設置為目標頻段的1.5倍再通過模式序號和場分布圖做篩選。具體參數(shù)方面我通常設置“所需模式數(shù)”為8到12搜索范圍是[0, 2×f_max]。f_max是目標頻率上限。如果范圍太窄高頻率模式會被漏掉如果范圍太寬會混入無效的數(shù)值模式。求解完成后COMSOL會在日志中給出“拒收特征值”列表這些被剔除的模式往往暗示著數(shù)值偽?；蜻吔缭O置問題值得仔細查看。3. 一維案例復現(xiàn)透射譜與一維禁帶3.1 一維多層膜模型的幾何參數(shù)設定一維光子晶體最常見的形式是交疊膜堆。我復現(xiàn)的一個典型算例是紫外波段的多層膜SiO2層與TiO2層交替排列厚度分別為95 nm和65 nm周期數(shù)10。在COMSOL里我選擇用二維模型來搭建幾何雖然結構是一維周期但二維模型能直觀觀察場分布也為后續(xù)斜入射計算留了余地。幾何構建時我會先用一個矩形代表整個膜堆再用“分割面”功能按層厚切成一系列子域。如果每一層都建獨立矩形再拼接后期改厚度會非常痛苦。分割面的操作在5.6里支持參數(shù)控制把層厚定義成全局參數(shù)后改一個數(shù)值整個幾何自動更新。材料方面SiO2折射率設為1.46TiO2設為2.35特別注意要把材料屬性里的損耗虛部清零否則禁帶位置會偏移。3.2 端口、周期邊界與入射波設置一維膜堆的透射和反射譜需要在結構兩側設置端口邊界條件。COMSOL的“端口”特性支持多模式設置入射端口的模式類型要選“衍射級”端口寬度必須包含一個完整周期。如果端口寬度小于一個周期透射率曲線會出現(xiàn)莫名其妙的震蕩這個問題非常隱蔽我調(diào)試了整整一天才找到原因。上下兩側需要設置Floquet周期性邊界把x方向的周期落實到模型中。注意上下邊界不能使用默認的PEC或者PBC否則會引入非物理的反射。頻率掃描范圍設置為320 nm到440 nm波長跑完結果后能看到反射譜在390 nm附近出現(xiàn)明顯的禁帶這是兩材料界面布拉格反射最強烈的波長位置。3.3 結果對照與網(wǎng)格精度控制我最初的版本誤差很大禁帶邊緣頻率比書中值偏移了7%左右。排查后發(fā)現(xiàn)是網(wǎng)格太粗65 nm厚的薄層里默認網(wǎng)格只剖了一層單元邊界處電磁場分布根本沒解析出來。把網(wǎng)格最大單元尺寸調(diào)整為25 nm后偏差縮小到了0.8%這個精度已經(jīng)滿足大多數(shù)工程需求。這個坑讓我養(yǎng)成了一個習慣所有薄膜結構每一層材料至少要跨4層網(wǎng)格。具體做法是使用“邊界層網(wǎng)格”在每層材料介質(zhì)界面處強制加密。多層結構用邊界層網(wǎng)格增加的自由度很少但對能帶位置的影響非常明顯?？梢哉f一維光子晶體仿真精度不夠大概率是網(wǎng)格的問題而不是求解器或物理設置的問題。4. 二維案例復現(xiàn)能帶結構中的TM/TE模式4.1 正方晶格與三角晶格的建模差異二維案例是這套案例庫中數(shù)量最多的部分。二維結構既能展現(xiàn)周期結構的共性又能通過不同的晶格排列得到豐富的能帶性質(zhì)。正方晶格和三角晶格的差別不僅僅在幾何排布上更關鍵的是它們的倒格子形狀和高對稱點路徑完全不同。正方晶格的倒格子仍是正方布里淵區(qū)高對稱路徑為Γ-X-M-Γ三角晶格的倒格子是六角對稱高對稱路徑為Γ-M-K-Γ。在幾何搭建時正方晶格只需一個正方形晶胞x和y方向設兩個周期性邊界。三角晶格則必須用平行四邊形晶胞兩個基矢長度相等但夾角120°。這里有一個常見錯誤很多人直接在正方形外框里放一個圓形空氣孔來模擬三角晶格這等于改變了晶格對稱性算出的能帶并不屬于真正的三角晶格。我的標準做法是建一個平行四邊形晶胞使用全局參數(shù)定義頂點坐標比如a600 nm、r180 nm然后利用三角函數(shù)關系算出平行四邊形的斜邊頂點。Floquet邊界恰好映射兩個基矢方向這樣才保證計算模型的對稱性正確。4.2 k路徑掃描的參數(shù)化實現(xiàn)二維案例能帶計算中k路徑掃描是最容易出錯也最耗時的環(huán)節(jié)。我在全局參數(shù)里定義了三段掃描變量s1、s2、s3分別對應Γ-M、M-K、K-Γ三段路徑。每一段路徑用線性插值把掃描參數(shù)映射到k點坐標。比如Γ-M段s從0到1kx從0映射到0.5ky從0映射到0這里的坐標都是相對于倒格子基矢的無量綱坐標。為了把三段路徑拼接到一次研究中我會定義一個總的掃描變量s通過分段函數(shù)判斷s落在哪一段再切換對應的k坐標表達式。這個寫法看起來繁瑣但在COMSOL中可以用“階梯函數(shù)”或“if條件表達式”實現(xiàn)設置完成后整個能帶掃描是一次性跑完的。4.3 參數(shù)掃描與禁帶優(yōu)化二維結構最有價值的應用就是禁帶優(yōu)化。比如設計工作在通信波長1550 nm附近的空氣孔光子晶體晶格常數(shù)a和空氣孔半徑r是兩個最關鍵的自由度。書里通常給一組基準參數(shù)但實際設計時需要掃描r/a比值。我在案例庫中準備了兩類參數(shù)掃描模型。第一類是固定a、掃描r觀察TM模禁帶寬度變化。典型結果是r/a從0.2增大到0.35時TM模禁帶逐漸變寬超過0.4后禁帶反而開始收窄因為空氣孔之間的介質(zhì)墻太薄高次模開始出現(xiàn)。第二類是固定r/a、整體縮放a觀察歸一化禁帶位置的變化。這類模型用來驗證光子晶體的縮放定律歸一化頻率a/λ基本保持不變這是周期結構設計的理論基礎也是能帶圖與實驗對照的關鍵參照。參數(shù)掃描時內(nèi)存占用不小。我在32 GB內(nèi)存的機器上一個三角晶格案例單次求解約2分鐘20組掃描約40分鐘。如果網(wǎng)格超過10萬自由度建議用“輔助掃描”來代替“參數(shù)掃描”能大幅減少內(nèi)存壓力。我實測下來兩種方式的精度差異可以忽略。5. 三維案例復現(xiàn)從幾何搭建到資源調(diào)配5.1 木堆結構的幾何布爾與域設置三維光子晶體案例中木堆結構非常經(jīng)典。它由多層介電柱堆疊而成每層柱子方向旋轉(zhuǎn)90度四層構成一個周期。在COMSOL中搭建木堆結構最讓人頭疼的是幾何布爾運算后的材料域標記。多個柱子做布爾并集后COMSOL有時會把交疊區(qū)域識別成內(nèi)部邊界導致后續(xù)網(wǎng)格無法跨邊界傳播。我的處理方式是在布爾運算之前給每根柱子做“顯式選擇”布爾操作時選擇“保留被選中的實體”把結構分為若干子域再逐個賦予材料。這樣即使后續(xù)做參數(shù)掃描修改柱寬材料分配也不會被打亂。每個三維案例我都會記錄幾何構建順序因為COMSOL的布爾運算是記錄在模型樹里的順序錯了后續(xù)修改幾乎無法進行。5.2 網(wǎng)格策略與內(nèi)存平衡三維能帶計算對硬件的需求很高。木堆結構如果直接使用默認的四面體網(wǎng)格單個晶胞就需要大概80萬到100萬個自由度內(nèi)存占用逼近16 GB。我的經(jīng)驗是兩步走先用粗網(wǎng)格快速試算獲得能帶的大致位置和模式數(shù)量然后再用細化網(wǎng)格在目標頻段精確計算。COMSOL 5.6的“自適應網(wǎng)格細化”在三維模型中很有效。它能自動識別場梯度大的區(qū)域在不增加整體網(wǎng)格數(shù)量的前提下修正局部精度。但自適應細化會增加迭代次數(shù)和應用時間所以只適合在最終求解階段開啟試算階段要保持關閉。內(nèi)存方面三維模型建議至少16 GB內(nèi)存求解時開啟多核并行速度提升非常明顯。5.3 三維場分布的后處理技巧三維案例除了能帶曲線通常還需要輸出漂亮的場分布圖。光子晶體場圖能直觀展示光與結構的相互作用線缺陷波導模式場集中在缺陷周圍木堆結構的光子禁帶模式場分布在介電柱之間的空隙里。在COMSOL里輸出場圖關鍵點是選對切面和位置。我常用的方法是用一個“工作平面”切過結構中間層疊加“高度圖”顯示場強配合透明顯示介質(zhì)結構。顏色表選擇需要注意彩色映射在黑白打印時會失真論文投稿建議使用灰度或雙色漸變。如果結果是復數(shù)場我會分別畫實部和模值實部用于觀察相位拓撲模值用于觀察能量分布。這兩種圖配合起來才能判斷模式是傳播態(tài)還是局域態(tài)。6. 進階弱形式方程求解色散光子晶體能帶6.1 內(nèi)置求解器為什么不夠用前面幾節(jié)的內(nèi)容都是基于“電磁波頻域”接口的“特征頻率”研究。這個方法對線性、無頻散、各向同性的介質(zhì)完全夠用。但遇到兩類結構內(nèi)置求解器就不太好辦了第一類是色散材料比如金屬、等離子體材料、增益介質(zhì)它們的介電常數(shù)隨頻率變化第二類是各向異性材料比如磁性光子晶體本構關系是張量形式?！疤卣黝l率”研究的求解流程是固定頻率后解出空間場。它需要先給定一個明確的介電常數(shù)值再算對應頻率。如果介電常數(shù)本身就是頻率的函數(shù)這個循環(huán)就變成了需要自洽求解的方程內(nèi)置求解器難以處理。而弱形式方程可以做到這也是“基于COMSOL弱形式方程求解色散光子晶體能帶”這個方向的初衷。6.2 弱形式方程的具體實現(xiàn)步驟弱形式的核心是把微分方程兩邊乘以一個任意檢驗函數(shù)再對整個求解域積分。以二維光子晶體TM模式為例控制方程是? × (1/ε_r(r) · ? × E_z) (ω2/c2) · E_z寫成等效的弱積分表達式之后被積函數(shù)中包含兩個部分第一項涉及電場梯度和檢驗函數(shù)梯度的乘積第二項是電場與檢驗函數(shù)相乘再乘上頻率項。在COMSOL的“弱形式PDE”接口里我們可以直接把被積函數(shù)寫成“Weak Expression”-(1/ε_r) * (Ex*test(Ex) Ey*test(Ey)) (ω2/c2) * Ez*test(Ez)這里的ε_r可以是任意表達式。比如Drude色散模型可以寫成ε_r ε_inf - ωp2/(ω2 i*γ*ω)其中ωp是等離子體頻率γ是阻尼率在COMSOL里定義為全局變量即可。此時介電常數(shù)將隨頻率和波矢量變化迭代求解時會自動滿足色散關系這正是弱形式方法相比內(nèi)置求解器的核心優(yōu)勢。操作上我會把E_z定義為弱形式PDE的因變量在“弱表達式”中輸入上面的被積函數(shù)然后添加全局方程來約束k點和本征頻率之間的關系。這樣在掃描k路徑時能帶結果能直接反映材料的頻散特性。6.3 不收斂問題的調(diào)試經(jīng)驗弱形式求解最大的注意點是初值。弱形式方程本質(zhì)是非線性的需要一個接近真實解的初始猜測。如果從零初值出發(fā)求解器幾乎必然發(fā)散。我的做法是先用無頻散模型算出同一結構的能帶把某個特征頻率作為色散模型的初始猜測再逐步增加色散項的強度。另一個問題是模式簡并。當兩個模式在同一個k點頻率相同弱形式求解器可能同時找到兩個解并發(fā)生混淆。我通常給結構加入一個很小的人工擾動比如把某個空氣孔半徑從180 nm改成180.1 nm專門破缺對稱性。這樣算出的能帶會有輕微劈裂但模式是干凈的。確認完模式屬性后再把擾動歸零重新計算。這個技巧在三維光子晶體中同樣有效尤其是處理重頻模式時特別管用。7. 高頻踩坑點與排查速查表7.1 能帶鋸齒網(wǎng)格問題還是物理問題能帶圖上出現(xiàn)鋸齒狀的尖刺多數(shù)情況下是網(wǎng)格問題但偶爾也有物理來源。判斷方法很簡單把目標頻段附近的網(wǎng)格尺寸減半重新計算同一區(qū)域。如果尖刺消失說明是網(wǎng)格欠分辨如果尖刺還在那可能是模式簡并劈裂需要用擾動法或模式分解來確認物理性質(zhì)。在二維案例中鋸齒經(jīng)常出現(xiàn)在布里淵區(qū)邊界附近這是因為邊界處場分布劇烈集中局部網(wǎng)格密度不足會帶來偽頻移。我專門在高對稱點附近加“角點細化”因為三角晶格中這些區(qū)域場梯度最大。這個方法在三維案例中也一樣有效能減少不少返工時間。7.2 模式缺失與特征值搜索范圍特征值缺失是復現(xiàn)案例時非常頭疼的問題。發(fā)現(xiàn)某個高對稱點附近缺少一條本該存在的能帶第一反應不該是改幾何而是檢查特征值搜索范圍。COMSOL特征值研究設置中有兩個關鍵參數(shù)一個是“所需模式數(shù)”決定求解器返回多少個特征模式另一個是“搜索基準點”決定搜索中心的頻率位置。如果頻帶跨度大我會把搜索基準點設置到目標頻帶中間而搜索范圍設置為整個感興趣頻段的1.5倍。同時把所需模式數(shù)提高到8以上跑完后手動過濾不需要的模式。還有一點務必注意COMSOL默認特征頻率單位是rad/s如果輸入以Hz為單位的值需要先做換算再填入。這個單位坑我踩過不止一次。7.3 歸一化頻率換算與單位陷阱光子晶體論文普遍用歸一化頻率a/λ而COMSOL的特征頻率輸出是freq單位rad/s。從freq換算到a/λ的公式很簡單a/λ a * freq / (2πc)其中c是真空光速。如果你習慣用波長作為輸入還要額外注意波長與頻域求解頻率之間的對應關系。這個換算看似簡單但做40多個案例時反復手動換算極容易出錯。我的解決方案是在“結果”節(jié)點中新建一個“全局變量探針”把換算公式直接定義成變量名作為一維繪圖的橫坐標。所有mph文件都采用這種輸出方式打開模型就能直接看到歸一化頻率的能帶圖免去了每次換算的工作量。7.4 常見問題速查表錯誤現(xiàn)象可能原因解決方向能帶曲線鋸齒形抖動網(wǎng)格過疏或界面網(wǎng)格不均勻加密界面網(wǎng)格開啟局部自適應細化禁帶位置偏移明顯材料折射率虛部未清零檢查材料屬性刪除損耗虛部特征值缺失頻率搜索范圍太窄增大搜索區(qū)間并調(diào)整所需模式數(shù)模式重復或發(fā)生交叉高對稱點簡并或?qū)ΨQ性太高加入人工微擾破缺對稱性結果全部為0端口或邊界條件錯誤檢查是否使用Floquet周期邊界而非通用周期邊界三維模型內(nèi)存溢出網(wǎng)格自由度過高關閉自適應細化簡化幾何或用輔助掃描弱形式方法不收斂初始猜測偏差過大先用無頻散能帶的計算結果做初值能帶圖整體變形k點坐標未做倒基矢變換重新設置Floquet邊界波矢量表達式透射率曲線震蕩端口寬度不等于一個完整周期將端口寬度調(diào)整為一個周期長度最后再分享一個小技巧mph文件保存前最好執(zhí)行一次“文件-壓縮模型”把求解過程中遺留的臨時數(shù)據(jù)和無效網(wǎng)格節(jié)點清理掉文件體積能縮小兩到三成。這不僅方便版本管理也方便和同行交換模型時減少傳輸壓力。另外復現(xiàn)案例時建議把每個模型的物理單位、頻段、材料參數(shù)記錄下來寫在README里否則三個月后回看你可能連自己的建模思路都忘了。這套案例庫目前仍在擴充下一步我準備加入更多拓撲光子晶體和谷態(tài)輸運的內(nèi)容有新的進展會繼續(xù)整理出來跟大家交流。