近場散射問題數(shù)值計算的超弱變分方法

欒 天，劉明輝

（1.北華大學 數(shù)學與統(tǒng)計學院，吉林 吉林132033；2.吉林大學 數(shù)學研究所，長春130012）

隨著近場光學理論的發(fā)展，近場散射問題受到人們廣泛關(guān)注.由于在遠場光學中存在分辨率衍射極限，而要突破這個極限獲取超分辨率的成像信息與探測信息，只能通過近場區(qū)域中的倏逝場實現(xiàn).因此，研究近場散射問題獲取倏逝波的信息十分重要.近場理論廣泛應用于納米技術(shù)、近場光學顯微鏡的研發(fā)和生物微樣本的無損成像等領(lǐng)域［1］.目前，關(guān)于近場時諧散射問題的數(shù)值研究結(jié)果報道較少，數(shù)值模擬多采用經(jīng)典的有限元方法、有限差分方法［2］和邊界積分方法［3］等.但隨著波數(shù)的增加，這些方法已不再適用.由于時諧波問題的解是振蕩的，因此當波數(shù)較高時，需要在單位波長距離內(nèi)選取足夠多的點才能得到相對穩(wěn)定的解.而且相位誤差的積累還會導致數(shù)值污染［4］.于是要達到給定的精度，就會迫使在每個波長內(nèi)使用更多的網(wǎng)格節(jié)點，計算量較大.基于此，本文采用另一類方法——Trefftz方法，即在區(qū)域剖分前提下，先在小單元上使用偏微分方程的解構(gòu)造離散空間，并在單元交界處借助數(shù)值方法使得離散解近似地滿足連續(xù)性；再將局部解“拼”接構(gòu)成一個整體近似解.這類方法主要包括單位分解法［5］、間斷強化法［6］、最小二乘法［7－8］和超弱變分方法［9］等.其中最典型的為超弱變分法.該方法源于區(qū)域分解技術(shù)，利用Green公式將問題轉(zhuǎn)化到網(wǎng)格邊界上處理.該方法目前已被應用于有界障礙物的散射問題中［10］.

本文針對近場光學中全內(nèi)反射顯微鏡的散射模型，提出了相應的超弱變分方法.首先，通過對區(qū)域網(wǎng)格剖分引入耦合參數(shù)，利用Green公式得到了與原邊值問題等價的超弱變分問題；其次，選取齊次Helmholtz方程的局部解（平面波，倏逝波）構(gòu)造基底，生成試探函數(shù)空間進行離散求解；最后，通過數(shù)值實驗驗證算法的有效性.結(jié)果表明，該算法能有效地數(shù)值模擬近場散射問題，適用于大波數(shù)情形，所需計算量少，收斂速度快.

1 數(shù)學模型

當沒有樣本時，記空間中的場分布為uref，稱其為參考場［11］.當在分界面Γ0上放置折射率為nS的樣本S時，參考場uref與樣本相互影響產(chǎn)生散射場us，如圖1所示.記此時空間中的總場u為

圖1 幾何模型Fig.1 Geometry of the model

根據(jù)Maxwell電磁理論，在TM極化情形下，u滿足如下方程：

空間中的波數(shù)為k0n（x），其中折射率

實際數(shù)值模擬中，需要將無界區(qū)域截斷，為此需引入人工Lipschitz邊界Γ，其外法方向記為ν，將問題限制在有界區(qū)域Ω內(nèi)，即Γ＝?Ω.為了保證所選取的解符合物理要求，要求散射場us在Γ上滿足吸收邊界條件：

于是，由式（1）可知總場u在Γ上滿足

從而所考慮的近場散射問題在數(shù)學上可描述為：已知參考場uref，求全場u滿足

其中g(shù)＝（?ν－i k0n（x））uref.

2 超弱變分問題

由于區(qū)域Ω中包含3種介質(zhì)，因此為了使全場u在不同介質(zhì)間的交界面處滿足連續(xù)性條件，需引入定義在剖分網(wǎng)格邊界上耦合參數(shù)：

由H的定義及u滿足邊值問題，有

將式（6）中兩式相減，得

將式（7）代入式（4），得

由u及其法向?qū)?shù)在Γk，j處的連續(xù)性及在Γ上滿足邊界條件，有

從而

將式（10）代入式（8）并關(guān)于k求和，得

引入算子F＝（Fk）∈L（X）定義為

其中ek是邊值問題

的解.

方程（12）稱為邊值問題（3）的超弱變分公式.

進一步，定義

按上述相同過程可推得式（11），與式（12）相減得

從而可推得式（9）成立，即u與其法向?qū)?shù)在內(nèi)邊界連續(xù)且滿足外邊界條件，于是由定義知，u是邊值問題（3）的解.

綜上，可得：

3 離散問題

選取有限維離散空間Xh?X代替超弱變分問題（14）中的空間X得到離散問題：求∈Xh，滿足

定義離散空間為

實際數(shù)值模擬中，需在每個小單元上選取方向互異的p個平面波函數(shù)和2個倏逝波函數(shù)作為ek，l.即若xk＝（xk，yk）表示單元Ωk的點，則

為簡單，平面波的方向dk，l通常選取如下：

由于在全反射過程中將出現(xiàn)倏逝波迅速衰減，因此為了能夠通過數(shù)值方法捕捉其信息，同時也選用如下倏逝波函數(shù)擴充平面波函數(shù)構(gòu)造基底：

選用上述兩個倏逝波，是基于本文的模型問題，依據(jù)參考場［11］uref選取的.實際應用中，可由具體模型考慮，選擇適合的倏逝波.此外，由于倏逝波是向上指數(shù)衰減的，因此在實際應用中并不要求在每個單元上都應用倏逝波函數(shù)構(gòu)造基底，而是僅在基座界面上方部分單元內(nèi)使用即可.

4 數(shù)值模擬

程序代碼使用 MATLAB語言.計算區(qū)域為Ω＝［－l，l］×［－l，l］，將Ω等矩形剖分成2nΩ×2nΩ（nΩ∈?）個小單元，單元邊長為l／nΩ.

例1 設k0＝1，l＝0.5，n＋＝45，n－＝75，nS＝60，θ＝π／4，則上半空間、下半空間和樣本介質(zhì)對應的波數(shù)分別為k＋＝45，k－＝75，kS＝60.

首先，數(shù)值求解參考場urefh.取定nΩ＝6，在每個單元上取p＝35個平面波函數(shù).由于倏逝波向上指數(shù)衰減，因此只在界面上方第一層單元內(nèi)引入倏逝波函數(shù)擴充平面波函數(shù)構(gòu)造基底.數(shù)值計算得到的參考場urefh實部等高線如圖2所示.

其次，為了驗證算法的有效性，考察算法關(guān)于p的收斂性.由于此時存在真解［11］：

從而可以計算真解uref和數(shù)值解在Ω上的相對誤差：選取初值p＝10，取nΩ＝3.當增加平面波數(shù)目，倏逝波數(shù)目保持不變，即增大p時.數(shù)值結(jié)果如圖3和圖4所示.圖3是相對誤差RE關(guān)于p的圖像，圖4是－lg（RE）關(guān)于p的圖像.由圖3和圖4可見，隨著基底數(shù)目的增加，誤差快速衰減，表明算法是關(guān)于p收斂的，而且收斂速度較快.當p＝50時，相對誤差已經(jīng)不超過1%.

最后，數(shù)值計算全場uh.取定nΩ＝6，p＝50.樣本大小為界面上方居中的兩個單元格.間接計算散射場＝uh－，其實部等高線如圖5所示.

圖2 urefh的實部等高線Fig.2 Contour for real part of urefh

圖3 相對誤差關(guān)于p的曲線Fig.3 Plot of relative error againstp

圖4 －lg（RE）關(guān)于p的曲線Fig.4 Plot of－lg（RE）againstp

圖5 ush的實部等高線Fig.5 Contour for real part of ush

綜上可見，本文算法在實際數(shù)值模擬中所需剖分單元少，收斂速度快，程序運行僅需數(shù)十秒，因此能有效計算近場散射問題.

［1］Carney P，Schotland J.Inside Out：Inverse Problems and Applications［M］.London：Cambridge University Press，2003.

［2］Farakawa H，Kawata S.Analysis of Image Formation in a Near－Field Scanning Optical Microscope：Effects of Multiple Scattering［J］.Opt Commun，1996，132（1／2）：170－178.

［3］Tanaka K，Tanaka M，Omoya T.Boundary Integral Equations for a Two－Dimensional Simulator of a Photon Scanning Tunneling Microscope［J］.J Opt Soc Amer A，1998，15（7）：1918－1931.

［4］Deraemaeker A，Babu?ka I，Bouillard P.Dispersion and Pollution of the FEM Solution for the Helmholtz Equation in One，Two and Three Dimensions［J］.J Numer Meth Eng，1999，46（4）：471－499.

［5］Babu?ka I，Melenk J M.The Partition of Unity Method［J］.Internat J Numer Meth Eng，1997，40（4）：727－758.

［6］Farhat C，Harari I，F(xiàn)ranca L P.The Discontinuous Enrichment Method［J］.Comput Meth Appl Mech Eng，2001，190（48）：6455－6479.

［7］Barnett A H，Betcke T.An Exponentially Convergent Nonpolynomial Finite Element Method for Time－Harmonic Scattering from Plygons［J］.SIAM J Sci Comput，2010，32（3）：1417－1441.

［8］LUAN Tian，WANG Yu－jie，ZHENG En－xi.Least Squares Method for Near－Field Scattering Problem ［J］.Journal of Jilin University：Science Edition，2013，51（5）：811－814.（欒天，王玉潔，鄭恩希.求解近場散射問題的最小二乘方法 ［J］.吉林大學學報：理學版，2013，51（5）：811－814.）

［9］Cessenat O，Després B.Application of an Ultra Weak Variational Formulation of Elleptic PDEs to Two－Dimensional Helmholtz Problem ［J］.SIAM J Numer Anal，1998，35（1）：255－299.

［10］Cessenat O，Després B.Using Plane Waves as Base Functions for Solving Time Harmonic Equations with the Ultra Weak Variational Formulation［J］.J Comput Acoust，2003，11（2）：227－238.

［11］Li P J.Numerical Simulations of Global Approach for Photon Scanning Tunneling Microscopy：Coupling of Finite－Element and Boundary Integral Methods［J］.J Opt Soc Amer A，2008，25（8）：1929－1936.