求解混合型擬周期散射問題的間斷Galerkin算法

欒 天, 王春艷, 郭 麗

(北華大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 吉林 吉林 132013)

1 擬周期散射問題

周期結(jié)構(gòu)的衍射問題在微光學(xué)領(lǐng)域應(yīng)用廣泛, 如納米級(jí)光學(xué)元件的設(shè)計(jì)、 光譜分析、 光纖通訊等[1]. 周期結(jié)構(gòu)通常也被稱為光柵, 關(guān)于光柵電磁學(xué)散射問題的數(shù)值研究目前已取得了很多成果[2-11]. 例如： Bao等[2]針對(duì)可穿透光柵衍射問題, 給出了一種自適應(yīng)完美匹配層(PML)有限元方法； Chen等[3]分析了PML算法在橫向電磁極化情形下的后驗(yàn)誤差估計(jì)； 尹偉石等[4]利用擬周期基本解, 應(yīng)用積分方程方法求解周期結(jié)構(gòu)表面散射問題, 并研究了該問題解的存在唯一性； Karimi等[5]考慮三維周期結(jié)構(gòu)的聲波外散射問題, 提出了一種有效的邊界元方法； 欒天等[6-7]采用最小二乘算法計(jì)算了可穿透的光柵衍射問題. 本文研究散射體為介質(zhì)和障礙混合型的擬周期散射問題, 其數(shù)學(xué)模型是關(guān)于Helmholtz方程的邊值問題.

圖1 散射體示意圖Fig.1 Schematic diagram of scatter

給定入射平面波ui, 產(chǎn)生散射場(chǎng)us, 若記u=ui+us為全場(chǎng), 則u滿足：

(1)

這里ki為散射體Ωi的波數(shù),k0為自由空間波數(shù). 邊界算子B可對(duì)應(yīng)Dirichlet,Neumann或Robin邊界條件.

u+-u-=?nu+-?nu-=0,

(2)

為保證數(shù)學(xué)問題的解存在唯一且符合物理要求, 還需散射場(chǎng)us滿足一定的輻射條件, 即散射場(chǎng)在無窮遠(yuǎn)處僅由有限多個(gè)向外傳播的平面波構(gòu)成. 設(shè)y0>0為充分大的實(shí)數(shù), 滿足矩形區(qū)域R=(-d/2,d/2)×(-y0,y0)包含D, 即D?R. 則在UR內(nèi)散射場(chǎng)us可由Rayleigh-Bloch展開[1], 即

(3)

于是, 本文所考慮散射問題的數(shù)學(xué)模型為: 給定入射平面波ui, 求全場(chǎng)u滿足式(1),(2), 且散射場(chǎng)us滿足輻射條件(3).

2 數(shù)值算法

算法區(qū)域剖分及符號(hào)如圖2所示. 首先, 將單元帶狀區(qū)域U剖分, 記

圖2 區(qū)域剖分及符號(hào)Fig.2 Domain decomposition and symbols

Dr+1=(-d/2,d/2)×(y0,+∞),

Dr+2=(-d/2,d/2)×(-y0,-∞).

其次, 在有界區(qū)域Di(i=1,2,…,q)內(nèi), 選取如下基本解函數(shù)的線性組合逼近散射場(chǎng):

(4)

(5)

同理, 對(duì)于區(qū)域D0, 選取區(qū)域Di(1,2,…,r)內(nèi)某閉曲線Ci上的點(diǎn)yi,j=(xi,j,yi,j)(i=1,2,…,r,j=1,2,…,Mi,Mi∈)作為點(diǎn)源, 定義逼近空間為

(6)

最后, 在Dr+1和Dr+2上采用Rayleigh-Bloch展開的有限項(xiàng)截?cái)嘧鳛閡s的近似, 這種近似可自然滿足有界外行平面波條件和擬周期條件. 在區(qū)域Dr+1上,

(7)

則相應(yīng)的近似空間為

同理, 對(duì)區(qū)域Dr+2, 有

綜上, 離散函數(shù)空間V定義為

由于選取的基函數(shù)都滿足Helmholtz方程, 因此可定義誤差匹配泛函為

其中[f]表示函數(shù)f在子區(qū)域公共邊界處的躍度, 定義為

從而算法的數(shù)值解uN為

(9)

利用對(duì)偶技巧可得如下誤差估計(jì):

‖u-uN‖0,Ω≤CJ(uN)1/2,

這里‖·‖0,Ω表示Ω上的L2范數(shù),J(uN)1/2=(J(uN))1/2.

定理1與文獻(xiàn)[12]中定理2.1的證明過程類似, 故略. 文獻(xiàn)[12]研究了折線形的良導(dǎo)體光柵, 本文所考慮的問題是具有擬周期結(jié)構(gòu)的介質(zhì)和障礙混合的散射體, 并且目標(biāo)泛函的定義與試探空間選取也不同. 定理1表明, 對(duì)任意非共振波數(shù)k0,J(uN)1/2控制解的內(nèi)部誤差, 因此可用于判斷算法的收斂性.

3 數(shù)值實(shí)驗(yàn)

下面利用MATLAB軟件對(duì)算法進(jìn)行數(shù)值模擬, 以驗(yàn)證算法在計(jì)算混合型擬周期散射問題時(shí)的有效性.

例1若單周期內(nèi)介質(zhì)散射體形狀為六葉草, 則其邊界曲線?D1參數(shù)方程為

Γ1(θ)=([0.8+0.1cos(6θ)]cosθ,[0.8+0.1cos(6θ)]sinθ);

若障礙散射體為卵形, 則其邊界曲線?D2的參數(shù)方程為

Γ2(θ)=([0.4+0.1cos(2θ)]cosθ,-2+[0.4+0.1cos(2θ)]sinθ).

其中0≤θ≤2π. 在Γ2上u=0, 即B為第一類邊界算子.

首先, 取y0=3將單周期帶狀區(qū)域U剖分, 如圖2所示. 然后選取點(diǎn)源構(gòu)造離散空間. 在D0內(nèi)選取點(diǎn)源

yj(θj)=([1+0.1cos(6θj)]cosθj,[1+0.1cos(6θj)]sinθj);

在D1,D2內(nèi)分別選取點(diǎn)

y1,j(θj)=([0.6+0.1cos(6θj)]cosθj,[0.6+0.1cos(6θj)]sinθj),

y2,j(θj)=([0.25+0.1cos(2θj)]cosθj,-2+[0.25+0.1cos(2θj)]sinθj).

為簡單, 取N1=M1=M2=p,p∈, 于是分別如圖2中綠色和紅色離散點(diǎn)所示.

其次, 選取入射場(chǎng)為平面波ui=exp{iαx-iβy}, 其中α=k0sinφ,β=k0cosφ, 入射角φ=-π/2, 自由空間波數(shù)k0=10, 介質(zhì)散射體內(nèi)波數(shù)k1=20, 周期d=3,p=30. 計(jì)算入射場(chǎng)ui、 散射場(chǎng)us和全場(chǎng)u, 其實(shí)部分別如圖3～圖5所示.

圖3 入射場(chǎng)實(shí)部Fig.3 Real part of incident field

圖4 散射場(chǎng)實(shí)部Fig.4 Real part of scattering field

最后, 分析p對(duì)匹配誤差J(uN)1/2的影響. 選取初始值p=10及步長為5, 逐漸增加p值, 計(jì)算相應(yīng)的泛函值J(uN)1/2, 數(shù)值結(jié)果如圖6所示. 由圖6可見, 隨著p的增加,J(uN)1/2快速衰減. 當(dāng)p=45時(shí), 精度已達(dá)到10-4, 表明收斂速度較快.

圖5 全場(chǎng)實(shí)部Fig.5 Real part of total field

圖6 J(uN)1/2收斂結(jié)果Fig.6 Convergence results of J(uN)1/2

數(shù)值算例結(jié)果表明, 該算法無需網(wǎng)格加密和空間截?cái)? 通過增加基底數(shù)目即能有效地計(jì)算全空間的散射場(chǎng), 從而檢驗(yàn)了本文算法在計(jì)算擬周期透射問題時(shí)的有效性.