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

方 宏 遠(yuǎn)， 林 皋＊， 張 蓓

（1.大連理工大學(xué) 建設(shè)工程學(xué)部 土木工程學(xué)院，遼寧 大連 116024；2.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；3.鄭州大學(xué) 水利與環(huán)境學(xué)院，河南 鄭州 450002）

0 引 言

近年來，電磁無損檢測技術(shù)被廣泛應(yīng)用于各個工程領(lǐng)域.很多工程結(jié)構(gòu)，比如道路、隧道等，都是層狀體系，利用電磁波對這些工程結(jié)構(gòu)進(jìn)行無損探測都可以歸結(jié)為電磁波在層狀體系中的傳播問題.分析電磁波在層狀介質(zhì)中的反射與透射，常采用傳遞矩陣方法[1]，但大部分工程材料都屬于有耗介質(zhì)，采用該方法計(jì)算有耗層狀體系，如果材料電導(dǎo)率比較大或者結(jié)構(gòu)層數(shù)比較多，容易出現(xiàn)數(shù)值失穩(wěn)問題，導(dǎo)致計(jì)算結(jié)果發(fā)散.

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

1 控制方程的推導(dǎo)[1，3]

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

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

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

式中：下標(biāo)s代表場量的水平分量x、y，下標(biāo)z表示場量的垂直分量.z＾表示z方向的單位矢量.則相應(yīng)的本構(gòu)參數(shù)變?yōu)?/p>

式中：ε′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 基于兩端邊值問題的精細(xì)積分方法[4－7]

精細(xì)積分是求解一階常微分方程的高精度算法，已被應(yīng)用于結(jié)構(gòu)動力分析、控制論、波導(dǎo)等諸多領(lǐng)域.該方法通過將原有積分步長再細(xì)分為2N個等長的微段，對于微段的層間矩陣?yán)脙缂墧?shù)展開求解.可以證明級數(shù)展開的截?cái)嗾`差低于計(jì)算機(jī)的舍入誤差，從而得到幾乎是計(jì)算機(jī)字長上精確的數(shù)值解.并利用結(jié)構(gòu)力學(xué)中的多層子結(jié)構(gòu)算法對這些微段進(jìn)行合并消元，在保證數(shù)值結(jié)果精度的前提下，大大提高了運(yùn)算效率.

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

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

將方程（10）對zb求導(dǎo)可得

在z＝zb處狀態(tài)方程（9）可以寫為

聯(lián)立方程（10）～ （12）可得

其中ea、hb可看作獨(dú)立無關(guān)的任意向量，所以有

當(dāng)za趨近于zb時，有邊界條件

2.1 區(qū)段方程的合并

對于相鄰兩區(qū)段[za，zb]、[zb，zc]，分別應(yīng)用方程（10）可以得到

對于合并后的區(qū)段[za，zc]，應(yīng)用方程（10）可以得到

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

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

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

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

2.2 區(qū)段矩陣E、F、G、Q的計(jì)算（2 N 類算法）

對于任意一層i，厚度為hi，將其劃分為等長的2N個小區(qū)間，N可以取得很大，例如20.由于間隔非常小，區(qū)段矩陣E、F、G、Q可按冪級數(shù)展開加以表示，并令τ＝hi/2N.當(dāng)τ足夠小時，冪級數(shù)取4項(xiàng)，截?cái)嗾`差已經(jīng)很小了.

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

將方程（23）代入方程（14），對比左右端項(xiàng)可得

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

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

3 邊界條件的推導(dǎo)

3.1 上部邊界條件

層狀體系如圖1所示，第1層為上部半無限空間，根據(jù)電磁場理論，第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的負(fù)特征值，代表上半空間中向下傳播的波.au＋＝ （au1au2）、au－＝（au3au4），為特征值相對應(yīng)的特征向量；Au＋＝（Au1Au2）T、Au－＝ （Au3Au4）T，為 相應(yīng)幅值；V＝ （exeyhxhy）T，為一四元列矢量.

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

定義反射系數(shù)矩陣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－為下半空間狀態(tài)方程中矩陣H的負(fù)特征值，代表下半無限空間中向下傳播的波.ad－＝ （ad3ad4），Ad－＝ （Ad3Ad4），符號表示與上半空間相同.

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

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

其中0為2×2零矩陣.

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

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

4 數(shù)值算例

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

例2 將例1中第1層介質(zhì)的電導(dǎo)率σ1變?yōu)?.05S/m，其他參數(shù)均不變，則垂直入射時，第1層與第2層界面處透射系數(shù)如表2所示.

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

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

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

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

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

5 結(jié) 論

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

[1] 鄭宏興，葛德彪.廣義傳播矩陣法分析分層各向異性材料對電磁波的反射與透射[J].物理學(xué)報(bào)，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] 趙麗濱，張建宇，王壽梅.精細(xì)積分方法的穩(wěn)定性和精度分析[J].北京航空航天大學(xué)學(xué)報(bào)，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] 周永祖.非均勻介質(zhì)中的場與波[M].聶在平，柳清伙，譯.北京：電子工業(yè)出版社，1992.Chew Weng－cho.Waves and Fields in Inhomogeneous Media [M].NIE Zai－ping，LIU Qing－h(huán)uo，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－h(huán)ao，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－h(huán)ao，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