有限元方法求解金屬襯底介質(zhì)片對(duì)平面波的反射問題

徐源，孫一月

(1.中國(guó)傳媒大學(xué) 理學(xué)院，北京100024；2.中國(guó)人民武裝警察部隊(duì)工程大學(xué) 理學(xué)院，陜西 西安710086)

有限元方法求解金屬襯底介質(zhì)片對(duì)平面波的反射問題

徐源1，孫一月2

(1.中國(guó)傳媒大學(xué) 理學(xué)院，北京100024；2.中國(guó)人民武裝警察部隊(duì)工程大學(xué) 理學(xué)院，陜西 西安710086)

任意平面波可以被分解為只有Z軸向電場(chǎng)的Ez-極化平面波和只有Z軸向磁場(chǎng)的Hz-極化平面波。我們需要研究的是金屬襯底介質(zhì)片對(duì)這兩種波的反射能力。根據(jù)標(biāo)量的赫爾姆霍茲方程分別得到相應(yīng)的極化方程，在方程的邊界上采用不同邊界條件。最后我們通過采用有限元的思想來(lái)分析問題，通過Matlab語(yǔ)言實(shí)現(xiàn)有限元方法的過程來(lái)解決這個(gè)問題。文章的最后給出了兩個(gè)數(shù)值實(shí)驗(yàn)來(lái)驗(yàn)證真解和數(shù)值解之間的誤差關(guān)系保證了算法的正確性和時(shí)效性。

極化平面波；有限元；Matlab

1 問題的描述

在圖1中，一束均勻的入射波在非均質(zhì)介質(zhì)片面上發(fā)生了反射。介質(zhì)片的寬度為L(zhǎng)，相對(duì)電導(dǎo)率為εr，相對(duì)的磁導(dǎo)率為μr都視為x的方程。周圍的自由空間介質(zhì)中εr=μr=1，接下來(lái)我們需要考慮的是介質(zhì)片的反射能力。

眾所周知任意平面波可以被分解為一個(gè)只有Z軸向電場(chǎng)的Ez-極化平面波和一個(gè)只有Z軸向磁場(chǎng)的Hz-極化平面波。因此我們可以充分考慮這兩種波的案例。對(duì)于Ez-極化平面波例子，其入射波的表達(dá)式可以是：

(1)

圖1 金屬介質(zhì)片平面波的反射圖

E0為一常量，指的是入射場(chǎng)的大小，θ代表圖1中的入射角度。顯而易見為了滿足在與x坐標(biāo)軸相垂直交界面處場(chǎng)的連續(xù)性條件，那么總場(chǎng)必須要有一個(gè)共同的系數(shù)e-jk0ysinθ。根據(jù)標(biāo)量的赫爾姆霍茲方程使得電場(chǎng)Ez滿足：

(2)

對(duì)于邊界條件我們采用的是齊次的狄利克雷邊界條件

Ez|x=0=0

(3)

相似地，對(duì)于極化波，入射波可以表示為：

(4)

總的磁場(chǎng)滿足標(biāo)量的赫爾姆霍茲方程：

(5)

對(duì)應(yīng)的紐曼邊界條件為：

(6)

2 解的分析

此問題的解決第一步需要把介質(zhì)片分成很多的薄層，定義為M層(圖2)每一層的相對(duì)電導(dǎo)率εr和相對(duì)磁導(dǎo)率可以近似地看作為常值，用εrm和μrm其中(m=1，2，3，…，M)。因此關(guān)于方程2在第m層的解可以表示為：

Ezm=(Amejkxmx+Bme-jkxmx)e-jk0ysinθ

(7)

圖2

我們強(qiáng)制Ez和Hy在兩層之間的交界面處是連續(xù)的，指出第m層和第(m+1)層的時(shí)候，我們有：

(8)

這里的xm+1指的是交界面處的位置，Rm=Bm/Am，并且有：

(9)

由于邊界條件式3，一旦有A1+B1=0成立，那么R1=-1。根據(jù)式8我們又可以計(jì)算出R2，R3等等。最后我們就可以計(jì)算出我們感興趣的反射系數(shù)RM+1，我們用|RM+1|2來(lái)表示介質(zhì)片反射能力的百分比。相類似的過程可以直接用在磁場(chǎng)的Hz極化波上面：

(10)

由邊界條件式6，一旦有A1-B1=0成立，那么R1=1。根據(jù)式8我們又可以計(jì)算出R2，R3等等。我們就可以計(jì)算出我們感興趣的Hz極化波反射系數(shù)RM+1。

3 有限元的解

首先，讓我們用有限元的方法來(lái)解決這個(gè)問題。首先確定問題的求解區(qū)域。顯然問題的求解區(qū)域是半無(wú)窮的(0≤x≤∞)，遺憾是有限元方法不能直接的用于這種半有界的求解域。于是我們需要適當(dāng)?shù)目s小下求解域，我們可以引進(jìn)適當(dāng)?shù)娜斯み吔鐥l件。為了滿足解的有效性，人工的邊界和求解域直接應(yīng)該足夠的接近。所以針對(duì)此問題最可能的選擇是令x=L，這里的L指的是板的厚度。

接下來(lái)要做的是在此邊界上找一個(gè)適當(dāng)?shù)倪吔鐥l件。需要指出的是在板的外部，總場(chǎng)可以表示為入射場(chǎng)和反射場(chǎng)的疊加。因此對(duì)于Ez極化波我們有：

Ez(x)=E0ejk0xcosθ+RE0e-jk0xcosθxgt;L

(11)

這里的R代表反射系數(shù)以及公共系數(shù)e-jk0xcosθ已經(jīng)表示過了，對(duì)Ez進(jìn)行對(duì)x求偏導(dǎo)數(shù)，我們有：

(12)

也可以表示為：

(13)

或者

(14)

此條件適用于在μr=1處介質(zhì)片外表面的邊界上。如果希望將邊界放在介質(zhì)片的內(nèi)表面，那么可以用到式14的邊界條件需要指出的是：

這是針對(duì)Ez和Hz連續(xù)邊界條件的結(jié)果，因此我們不難發(fā)現(xiàn)：

(15)

作為邊界條件僅僅在介質(zhì)片的內(nèi)表面是可行的。顯然式13和式14是等價(jià)的，它們的選擇并不會(huì)影響方程的解。相類似地針對(duì)極化波我們同樣可以找到適當(dāng)?shù)倪吔鐥l件強(qiáng)加在Hz如：

(16)

或

(17)

對(duì)于式2和式5的不同，我們分別給出了在x=0處的邊界條件和在x=L處推導(dǎo)出的邊界條件。我們就可以繼續(xù)運(yùn)用有限元方法來(lái)求解問題的數(shù)值解。我們首先是把介質(zhì)片或者是區(qū)域(0，L)分成了一定數(shù)量的薄層，再次將M層每一層都看作是一個(gè)單元。因此我們就可以通過編寫有限元方法的計(jì)算機(jī)語(yǔ)言可以高效地解決這個(gè)問題。一旦場(chǎng)在x=L處被解決了，那么方程11的反射系數(shù)也可以得到，或者更具體地對(duì)于Ez極化波的反射系數(shù)來(lái)說：

Hz極化波來(lái)說也是類似，就是把Ez用Hz來(lái)代替，E0用H0來(lái)代替。

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

為了保證有限元解的有效性，我們令電導(dǎo)率εr=4+(2-0.1j)(1-x/L)2，相對(duì)應(yīng)的磁導(dǎo)率μr=2-0.1j。對(duì)于介質(zhì)片的厚度采取五倍于自由空間的波長(zhǎng)。實(shí)驗(yàn)的結(jié)果已經(jīng)展示在圖3。采用50單元(51個(gè)節(jié)點(diǎn))和100單元(101個(gè)節(jié)點(diǎn))的有限元的數(shù)值解和真解的比較。正如我們所看到的節(jié)點(diǎn)的數(shù)量越多，數(shù)值解就越接近真解。

5 實(shí)驗(yàn)分析

通過Matlab語(yǔ)言實(shí)現(xiàn)的有限元算法來(lái)求解金屬襯底介質(zhì)板對(duì)平面波的反射問題，我們不難發(fā)現(xiàn)，隨著求解單元數(shù)量的增加，數(shù)值解與真實(shí)解之間的誤差會(huì)越來(lái)越小。這達(dá)到了我們預(yù)先的估計(jì)和目的，用有限元的方法來(lái)解決偏微分方程問題的時(shí)效性和準(zhǔn)確性優(yōu)勢(shì)明顯。根據(jù)完整的有限元的求解過程構(gòu)編寫程序語(yǔ)言是本文的最終目的。

(a)電場(chǎng)極化算例 (b)磁場(chǎng)極化算例圖3 金屬襯底介質(zhì)片的反射系數(shù)εr=4+(2-0.1j)(1-x/L)2，μr=2-0.1j，L=5λ

[1]Harris，Edward G.Introduction to modern theoretical physics[M].V I，1975.

[2]Jianming Jin.The finite element method in electromagnetics[M].Canada：John Wiley amp; Sonsinc，2002：83-90.

[3]嚴(yán)登俊，黃學(xué)良，胡敏強(qiáng)，勾磊.有限元網(wǎng)格生成技術(shù)分析[J].微特電機(jī)，1999(01).

[4]傅為農(nóng)，江建中.關(guān)于電磁場(chǎng)數(shù)值分析的若干認(rèn)識(shí)[J].微特電機(jī)，1998(05).

[5]劉春太，楊曉東，申長(zhǎng)雨，陳靜波.任意平面區(qū)域的變尺寸有限元網(wǎng)格劃分[J].計(jì)算力學(xué)學(xué)報(bào)，2000(01).

[6]Gerald Recktenwald.數(shù)值方法和MATLAB實(shí)現(xiàn)與應(yīng)用[M].北京：機(jī)械工業(yè)出社，2004.

[7]Snyder A W，Love J D.光波導(dǎo)理論[M].北京：人民郵電出版社，1991.

(責(zé)任編輯：宋金寶)

FiniteElementMethodtoSolvePlane-waveReflectionbyaMetal-backedDielectricSlab

XU Yuan1，SUN Yi-yue2

(1.School of Science，Communication University of China，Beijing 100024；2.College of Science，Engineering University of the Chinese People’s Armed Police Force，Xi’an 710086)

It is well known that any plane wave can be decomposed into an Ez-polarized plane wave having only a z-compoment for the electric field and a Hz-polarized plane wave having only a z-compoment for the magnetic field，we need to find the power reflected by the slab for two polarization cases only.With observation，the scalar Helmholtz equation governing the functions，the boundary condition for equations is different.We consider the finite element solution to this problem and write Matlab code to solve the problem. Finally this dissertation gives the running result of this system and analyzes the accuracy of the algorithm and its stabilization.

polarization plane wave；finite element solution；Matlab

O24

A

1673-4793(2017)06-0032-04

2017-07-06

徐源(1990-)，男(漢族)，山東臨沂人，中國(guó)傳媒大學(xué)碩士研究生.E-mail：xuyuan@cuc.edu.cn