求解電磁波在層狀有耗介質中反射和透射的精細積分方法

方 宏 遠， 林 皋＊， 張 蓓

（1.大連理工大學 建設工程學部 土木工程學院，遼寧 大連 116024；2.大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024；3.鄭州大學 水利與環境學院，河南 鄭州 450002）

0 引 言

近年來，電磁無損檢測技術被廣泛應用于各個工程領域.很多工程結構，比如道路、隧道等，都是層狀體系，利用電磁波對這些工程結構進行無損探測都可以歸結為電磁波在層狀體系中的傳播問題.分析電磁波在層狀介質中的反射與透射，常采用傳遞矩陣方法[1]，但大部分工程材料都屬于有耗介質，采用該方法計算有耗層狀體系，如果材料電導率比較大或者結構層數比較多，容易出現數值失穩問題，導致計算結果發散.

精細積分方法[2]是一種高精度且無條件穩定[2]的數值算法，很適合求解層狀體系中電磁波的反射和透射問題.為了避免出現傳遞矩陣方法中指數矩陣相乘的問題，可采用基于兩端邊值問題的精細積分方法，將控制方程按邊值問題進行計算，從而避免了在計算過程中出現指數增大的情況，可以保持數值結果的精度和穩定性.

1 控制方程的推導[1，3]

不失一般性，文中介質按各向異性考慮.此時，頻域麥克斯韋方程組可以表述為

式中：h表示磁場向量；e表示電場向量；μ為磁導率；ε′為復介電常數，ε′＝ε＋σ/iω，其中ε是介電常數，σ是電導率.各向異性介質中μ、ε′均為二階張量.

由電磁場理論可知，在兩種介質界面處，電場和磁場的水平分量連續.為了便于在分層界面上匹配邊界條件，可將旋度算子和場量按下式進行分解：

式中：下標s代表場量的水平分量x、y，下標z表示場量的垂直分量.z＾表示z方向的單位矢量.則相應的本構參數變為

式中：ε′s、μs是2×2的矩陣，ε′sz、μsz是2×1的矩陣，ε′zs、μzs是1×2的矩陣.

將式（3）、（4）代入式（1）、（2）整理可得

電場和磁場的頻域列式可以寫為

將式（7）、（8）代入式（5）、（6），控制方程可以寫為

式中：es＝ （exey）T，hs＝ （hxhy）T，H11、H12、H21、H22均為2×2的矩陣.

2 基于兩端邊值問題的精細積分方法[4－7]

精細積分是求解一階常微分方程的高精度算法，已被應用于結構動力分析、控制論、波導等諸多領域.該方法通過將原有積分步長再細分為2N個等長的微段，對于微段的層間矩陣利用冪級數展開求解.可以證明級數展開的截斷誤差低于計算機的舍入誤差，從而得到幾乎是計算機字長上精確的數值解.并利用結構力學中的多層子結構算法對這些微段進行合并消元，在保證數值結果精度的前提下，大大提高了運算效率.

邊值問題精細積分方法并沒有采用傳統的差分格式對控制方程（9）進行離散求解.對于線性系統中任意區段[za，zb]，兩端處電場和磁場的關系一定可以表示為

式中：ea、eb、ha、hb分別表示兩端的電場和磁場向量.E、F、G、Q為待求的2×2矩陣，均為za、zb的函數.

將方程（10）對zb求導可得

在z＝zb處狀態方程（9）可以寫為

聯立方程（10）～ （12）可得

其中ea、hb可看作獨立無關的任意向量，所以有

當za趨近于zb時，有邊界條件

2.1 區段方程的合并

對于相鄰兩區段[za，zb]、[zb，zc]，分別應用方程（10）可以得到

對于合并后的區段[za，zc]，應用方程（10）可以得到

將方程（17b）代入方程（16a）求出

將方程（16a）代入方程（17b）得

將方程（19）、（20）代入式（16b）、（18a）整理可得

對比方程（18）和（21）可得

2.2 區段矩陣E、F、G、Q的計算（2 N 類算法）

對于任意一層i，厚度為hi，將其劃分為等長的2N個小區間，N可以取得很大，例如20.由于間隔非常小，區段矩陣E、F、G、Q可按冪級數展開加以表示，并令τ＝hi/2N.當τ足夠小時，冪級數取4項，截斷誤差已經很小了.

式中：θi、γi、φi、ψi（i＝1，2，3，4）均為2×2矩陣，I為2×2單位矩陣.

將方程（23）代入方程（14），對比左右端項可得

這里需要注意的是，τ特別小，導致F′（τ）、E′（τ）特別小，如果直接將其與單位陣相加的話，會導致計算精度喪失，所以不能直接利用式（22）進行組合，需改寫為以下形式：

由于同一種介質層內所有區段均相等，采用式（24）消元一次，區段數就會減少一半，經過N次消元后，即可求出這一層的區段矩陣.依此類推可以求得每一介質層的區段矩陣Ei、Fi、Gi、Qi.然后按照式（22）組裝整個層狀體系的區段矩陣Ec、Fc、Gc、Qc.

3 邊界條件的推導

3.1 上部邊界條件

層狀體系如圖1所示，第1層為上部半無限空間，根據電磁場理論，第1層中傳播的波可以表述為向上傳播的波和向下傳播的波兩部分.

式中：exp （iλu＋（z－z0））＝ diag｛exp （iλu1（zz0））， exp（iλu2（z－z0））｝，λu＋為方程（9）中矩陣H的正特征值，代表上半空間中向上傳播的波.exp（iλu－（z－z0））＝diag｛exp （iλu3（z－z0）），exp（iλu4（z－z0））｝，λu－為矩陣H的負特征值，代表上半空間中向下傳播的波.au＋＝ （au1au2）、au－＝（au3au4），為特征值相對應的特征向量；Au＋＝（Au1Au2）T、Au－＝ （Au3Au4）T，為 相應幅值；V＝ （exeyhxhy）T，為一四元列矢量.

圖1 層狀體系示意圖Fig.1 Sketch of layered system

定義反射系數矩陣R，R＝Au＋/Au－，R為2×2矩陣，式（25）可重寫為其中I為2×2單位矩陣.

在z＝z0處將上式展開，從而得到上部半空間的邊界條件

式中au11、au12、au21、au22均為2×2矩陣，可由au分塊獲得.

3.2 下部邊界條件

第n＋1層為下部半無限空間，下部半空間只有向下傳播的波，所以

式中：exp （iλd－（z－zn））＝ diag ｛exp （iλd3（zzn）），exp（iλd4（z－zn））｝，λd－為下半空間狀態方程中矩陣H的負特征值，代表下半無限空間中向下傳播的波.ad－＝ （ad3ad4），Ad－＝ （Ad3Ad4），符號表示與上半空間相同.

定義折射系數矩陣T，T＝Ad－/Au－，T為2×2矩陣.

所以下半空間的波可以寫為

其中0為2×2零矩陣.

在z＝zb處，展開上式可得下半空間的邊界條件

將式（27）、（30）代入方程（10），即可求得層狀體系的反射系數和折射系數.這里的E、F、G、Q是采用精細積分方法求得的整個層狀體系的區段矩陣Ec、Fc、Gc、Qc.

4 數值算例

iωε′，其余元素均為0.其中為電磁波在空氣中的傳播常數，θ為入射角.設第1層介質介電常數ε1＝9ε0，電導率σ1＝0.05S/m，厚度d1＝0.3m；第2層介質介電常數ε2＝30ε0，電導率σ2＝0.05S/m，為半無限空間.兩層介質磁導率均等于真空磁導率，則入射角為0°和30°時，表面反射系數如表1所示.

例2 將例1中第1層介質的電導率σ1變為4.05S/m，其他參數均不變，則垂直入射時，第1層與第2層界面處透射系數如表2所示.

例3 例1中其他參數不變，第1層的介電常數按各向異性考慮，εxx＝9ε0，εxy＝8ε0，εxz＝7ε0，9ε0，εzz＝10ε0，則空氣與第1層介質界面處反射系數如表3所示.

例1 均勻平面波由空氣入射一個兩層體系，設每層介質均為各向同性，空氣層的介電常數和磁導率分別為真空中的介電常數和磁導率ε0、μ0.那 么 方 程 （9）中 矩 陣H可 以 簡 化 為

表1 兩層體系反射系數對比結果Tab.1 Comparison result of reflection coefficient in two layered medium

表2 兩層體系透射系數對比結果Tab.2 Comparison result of transmission coefficient in two layered medium

表3 各向異性兩層體系反射系數Tab.3 Reflection coefficient of two layered anisotropic medium

5 結 論

對比精細積分方法和廣義傳遞矩陣法數值結果可知，對于低耗層狀體系，兩種方法都能獲得精確的結果，但如果介質電導率比較大，廣義傳遞矩陣方法的數值結果就出現了不穩定性.這主要是由于該方法求解結構層總傳遞矩陣過程是一個指數矩陣相乘的問題，會出現大數溢出的現象，電導率越大，傳遞矩陣中數值增大的項增長越快，越容易造成數值不穩定.而精細積分方法中所示的層間矩陣合并過程，是一個指數矩陣相除的過程，避免了出現大數溢出的情況，保持了該算法的數值穩定性.利用基于兩端邊值問題的精細積分方法計算層狀體系反射系數和透射系數，具有更高的精度和數值穩定性.

[1] 鄭宏興，葛德彪.廣義傳播矩陣法分析分層各向異性材料對電磁波的反射與透射[J].物理學報，2000，49（9）：1702－1706.ZHENG Hong－xing，GE De－biao.Electromagnetic wave reflection and transmission of anisotropic layered media by generalized propagation matrix method[J].Acta Physica Sinica，2000，49（9）：1702－1706.（in Chinese）

[2] 趙麗濱，張建宇，王壽梅.精細積分方法的穩定性和精度分析[J].北京航空航天大學學報，2000，26（5）：569－572.ZHAO Li－bin，ZHANG Jian－yu，WANG Shou－mei.Stability and precision analysis for precise integration method [J]. Journal of Beijing University of Aeronautics and Astronautics，2000，26（5）：569－572.（in Chinese）

[3] 周永祖.非均勻介質中的場與波[M].聶在平，柳清伙，譯.北京：電子工業出版社，1992.Chew Weng－cho.Waves and Fields in Inhomogeneous Media [M].NIE Zai－ping，LIU Qing－huo，tran.Beijing：Electronic Industry Press， 1992. （in Chinese）

[4] GAO Qiang，ZHONG Wan－xie，Howson W P.A precise method for solving wave propagation problems in layered anisotropic media [J]. Wave Motion，2004，40（3）：191－207.

[5] ZHONG Wan－xie，LIN Jia－hao，GAO Qiang.The precise computation for wave propagation in a stratified materials [J].International Journal for Numerical Methods in Engineering，2004，60（1）：11－25.

[6] GAO Qiang，LIN Jia－hao，ZHONG Wan－xie，etal.A precise numerical method for Rayleigh waves in a stratified half space [J].International Journal for Numerical Methods in Engineering，2006，67（6）：771－786.

[7] ZHONG Wan－xie.Combined method for the solution of asymmetric Riccati differential equations [J].Computer Methods in Applied Mechanics and Engineering，2001，191：93－102