本發(fā)明屬于計算電磁學技術領域，具體涉及一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法。

背景技術：

時域有限差分(finite-differencetime-domain，fdtd)方法已經成功地應用于求解磁化等離子體中的電磁問題。由于它是一種顯式的數(shù)值計算的方法，其時間步長受柯西穩(wěn)定性條件的限制，不能取的太大，方法時效性不好，計算速度較慢，尤其在磁化等離子體中的電磁計算中，時間步長不僅僅受柯西穩(wěn)定性條件的限制，而且還與等離子體參數(shù)和差分方案有關，如此更進一步的限制了時間步長的選取，其計算速度和計算效率較空氣中的計算更差。為了消除柯西穩(wěn)定性條件的限制，人們提出了無條件穩(wěn)定時域有限差分方法，比如：交替方向隱式(alternating-direction-implicit，adi)的時域有限差分(adi-fdtd)方法和基于加權拉蓋爾多項式的時域有限差分(finite-differencetime-domainwithweighted-laguerre-polynomials，wlp-fdtd)方法。在這些方法中，adi-fdtd方法在使用較大的時間步長時會產生很大的色散誤差，而wlp-fdtd方法既能消除柯西穩(wěn)定性條件的限制，又能解決adi-fdtd方法在使用較大的時間步長時會產生很大的色散誤差這個難題，因此wlp-fdtd方法可以高效的求解磁化等離子體中的電磁問題。然而，這種wlp-fdtd方法在求解電磁場過程中，會產生一個大型的稀疏矩陣方程，直接求解此方程會使得計算較復雜，內存消耗較大，于是提出了一種因式分裂的wlp-fdtd方法，該方法在計算速度和計算效率上得到了較大的提高，但是由于該方法是通過添加微擾項，進行因式分裂得來的，計算中會產生分裂誤差，為了減小分裂誤差、提高計算精度，同時保證較快的計算速度，提出了一種高精度迭代的wlp-fdtd方法。

另外由于計算機容量的限制，電磁場的計算只能在有限區(qū)域進行。為了能模擬開域電磁波傳播過程，必須在計算區(qū)域的截斷邊界處給出吸收邊界條件。有人提出了完全匹配層(perfectlymatchedlayer，pml)吸收邊界，后來pml被廣泛應用于計算區(qū)域的截斷，而且被證明是非常有效的，但是研究發(fā)現(xiàn)這種傳統(tǒng)pml對低頻以及凋落波的吸收效果并不理想；使用帶有復頻率偏移(complexfrequencyshift，cfs)因子的擴展坐標的pml吸收邊界可以有效地改善傳統(tǒng)pml對低頻，凋落波的吸收效果。

技術實現(xiàn)要素：

本發(fā)明的目的是提供一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，計算精度高、計算速度快，且對于低頻和凋落波具有很好的吸收效果。

本發(fā)明所采用的技術方案是，一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，具體按照以下步驟實施：

步驟1：輸入模型文件；

步驟2：初始化參數(shù)和設置參數(shù)；

步驟3：添加場源到x方向上的電場分量系數(shù)中，設置磁場分量系數(shù)和分別記為初始場值和其中，添加的場源的表達式為：

jx(t)＝exp(-(t-t0)²/τ²)

其中，t表示時間，單位為秒，t0,τ為場源參數(shù)；

步驟4：更新計算整個計算區(qū)域的x方向上電場分量系數(shù)

步驟5：更新計算整個計算區(qū)域的y方向上電場分量系數(shù)

步驟6：更新計算整個計算區(qū)域的磁場分量系數(shù)

步驟7：更新計算整個計算區(qū)域的磁場分量系數(shù)

步驟8：將k+1賦值給k，并判斷迭代次數(shù)k是否達到預設值，若未達到預設值，則返回步驟4，若達到預設值，則執(zhí)行步驟9；

步驟9：更新計算整個計算區(qū)域的極化電流密度系數(shù)

步驟10：更新計算整個計算區(qū)域的電磁場分量系數(shù)的輔助變量；

步驟11：更新計算觀測點處的電磁場分量；

步驟12：將q+1賦值給q，并判斷拉蓋爾多項式的階數(shù)q是否達到預設值，若未達到預設值，則返回步驟3，若達到預設值，則結束。

本發(fā)明的特點還在于：

步驟1具體為：

計算區(qū)域大小nz，nz為z軸方向上網格數(shù)；空間步長δz；時間步長δt；真空中的電導率σ，磁導率μ0，介電常數(shù)ε0；等離子體碰撞頻率υ，等離子體頻率ωp，電子回旋頻率ωb；等離子體在計算區(qū)域中的位置；吸收邊界層數(shù)npml與相關參數(shù)κzmax，αzmax，σzmax；κzmax取整數(shù)，κzmax取值范圍為[1，60]；αzmax取值范圍為[01)；σzmax/σopt取值范圍為(0，12]，σopt＝(m+1)/150πδz，m取值范圍為[1，20]，δz取值范圍為λ為源的波長；仿真計算時長tf；迭代次數(shù)k，k≥0且為整數(shù)；加權拉蓋爾多項式的階數(shù)q，q≥0且為整數(shù)；時間尺度因子s，s取值范圍為[10⁹，10¹³]；觀測點；場源參數(shù)。

步驟2初始化具體為：

將整個計算區(qū)域的電磁場分量系數(shù)整個計算區(qū)域的極化電流密度系數(shù)整個計算區(qū)域的電磁場分量系數(shù)的和整個計算區(qū)域的極化電流密度系數(shù)的和整個計算區(qū)域的輔助變量和其中表示輔助變量表示輔助變量拉蓋爾多項式其中全部初始化為零；磁化等離子體參數(shù)(p1,p2,p3)初始化為p1＝0，p2＝0，pml系數(shù)(c1z，c2z，c3)初始化為c1z＝2/ε0s，c2z＝c1z，c3＝ε0/μ0，其中，ε0是空氣中的介電常數(shù)，μ0是空氣中的磁導率；

設置的參數(shù)具體為：

設置帶有cfs因子的sc-pml吸收邊界的參數(shù)σz,κz,αz；具體為：

σz＝σzmax|z-z0|^m/d^m；

κz＝1+(κzmax-1)|z-z0|^m/d^m；

αz＝αzmax；

式中z0為pml層與非pml截面位置，d是pml吸收邊界的厚度，κzmax取整數(shù)，κzmax取值范圍為[1，60]；αzmax取值范圍為[0，1)；σzmax根據σopt來設置，σzmax/σopt取值范圍為(0，12]；σopt＝(m+1)/150πδz，m取值范圍為[1，20]，δz取值范圍為λ為源的波長；

設置pml系數(shù)c1z，c2z；具體為：

c1z＝1/(κzαz+σz+0.5κzε0s)，c2z＝(2αz/ε0s+1)c1z；

設置等離子體參數(shù)為：

步驟4具體為：

步驟4.1：電場分量系數(shù)在計算區(qū)域的方程為：

其中，k表示迭代次數(shù)，n表示z軸上第n個計算網格處的位置，n+1/2表示z軸上第n個半網格處的位置，c2z|n表示系數(shù)c2z在z軸方向上第n個網格處的值，c2z|n-1/2表示系數(shù)c2z在z軸方向上第n-1個半網格處的值，δz|n表示z軸方向上第n個網格的大小，第表示第k+1次迭代第n-1個網格處的電場分量系數(shù)的值，

步驟4.2：使用追趕法求解步驟4.1的方程，得到整個計算區(qū)域的電場分量系數(shù)

步驟5具體為：

步驟5.1：電場分量系數(shù)在計算區(qū)域的方程為：

其中，表示第k+1次迭代第n-1個網格處的電場分量系數(shù)的值，

步驟5.2：使用追趕法求解步驟5.1的方程，得到整個計算區(qū)域的y方向上電場分量系數(shù)

步驟6的更新方程具體為：

其中

步驟7的更新方程具體為：

其中

步驟9更新公式具體為：

步驟10更新公式具體為：

其中fζz代表exz,eyz,hxz,hyz，因此表示輔助變量表示電磁場分量系數(shù)的和對z的求導，表示輔助變量的和，表示輔助變量

步驟11的更新公式具體為：

上式中u表示電磁場分量ex,ey,hx,hy，r表示電磁場分量的位置，u^q表示q階電磁場分量系數(shù)，是q階加權拉蓋爾多項式，是帶有時間尺度因子s>0的擴展時間，是q階拉蓋爾多項式。

本發(fā)明的有益效果是：

1).本發(fā)明一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，在直角坐標系下，通過用加權拉蓋爾多項式表示電磁場分量，來解時域麥克斯韋方程，使得在更新計算整個計算區(qū)域的電磁場分量系數(shù)時不涉及到時間步長，只是在最后計算觀測點處的電磁場分量時用到時間步長，因此計算過程中時間步長可以取得比柯西穩(wěn)定性條件限制的時間步長更大；

2).本發(fā)明一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，在求解求解磁化等離子體中的電磁場分量系數(shù)時，使用迭代的因式分裂的方案將大型稀疏矩陣方程分裂成兩個迭代的三對角矩陣方程，使得它在計算時比wlp-fdtd方法更簡單、計算速度更快、內存消耗更少而且比無迭代的因式分裂的wlp-fdtd方法精度更高；

3).本發(fā)明一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，在設置pml系數(shù)時，由于采用了cfs因子，并且通過調整cfs因子中的參數(shù)，可以使得該吸收邊界對低頻與凋落波的吸收更加有效；

4).本發(fā)明一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，由于采用了復擴展坐標系，使得pml在實現(xiàn)時避免了場的分裂且與媒質無關。

附圖說明

圖1是本發(fā)明實現(xiàn)方法的流程示意圖；

圖2是使用本發(fā)明的方法與解析解、wlp-fdtd方法和無迭代的因式分裂的wlp-fdtd方法計算的左旋極化波的反射系數(shù)振幅圖；

圖3是使用本發(fā)明的方法與解析解、wlp-fdtd方法和無迭代的因式分裂的wlp-fdtd方法計算的左旋極化波的反射系數(shù)相位圖。

具體實施方式

下面結合附圖和具體實施方式對本發(fā)明進行詳細說明。

本發(fā)明的一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，原理為：首先導出一維磁化等離子體中，電磁場所滿足的復擴展坐標系下的麥克斯韋方程，然后使用一維磁化等離子體中迭代的因式分裂的wlp-fdtd方法推導出整個計算區(qū)域的電磁場分量系數(shù)更新方程，最后采用公式(21)求解觀測點處的電磁場分量。

一維磁化等離子體中，電磁場所滿足的復擴展坐標系下的麥克斯韋方程推導過程如下：

在各向異性色散介質碰撞磁化等離子體中，麥克斯韋方程組和相關的聯(lián)立方程為

式中，h是磁場強度；e是電場強度；j是極化電流密度；ε0、μ0分別為真空中的介電常數(shù)和磁導率；υ是等離子體碰撞頻率；是等離子體角頻率的平方；ωb＝eb0/m是電子回旋頻率；b0是外磁場；e、m分別是電子的電量和質量。

設外磁場的方向為+z軸，方程(3)可寫為

應用擴展坐標的cfs-pml，僅考慮二維tez的情況，上述麥克斯韋方程組和相關的聯(lián)立方程可化為：

式中是帶有cfs因子的坐標伸縮因子，取為：

下面我們引入幾個輔助變量如下：

將(12)式分別代入(13)-(16)，然后做如下變化于是可得下面幾式：

我們可以知道

式中u代表ex、ey、hx、hy，是加權拉蓋爾多項式，q≥0；t≥0是q階拉蓋爾多項式。

將(21)、(22)代入(17)-(20)，然后應用加勒金的測試過程，可得

式中

s>0是時間尺度因子，q是加權拉蓋爾多項式的階數(shù)。

將(21)、(22)代入(6)-(11)，再應用加勒金的測試過程得到：

式中，

將(36)代入(37)得到

將(37)代入(36)得到

將(39)代入(32)得到

將(38)代入(33)得到

將方程(40)、(41)、(34)和(35)放置一起

整理后得到

其中

將(43)式寫成矩陣的形式

(i-a-b)x^q＝v^q-1(45)

式中

b＝p1p2/(1+p2),a1＝c2zdz/(1+p2),a2＝c2zc3dz

添加微擾項到(45)式得到

為了強調迭代的特點，將(47)式重寫為

引進中間變量于是(48)式可以分裂為

將上式展開得到

上兩式解得

將展開得到

將(52)、(53)和(54)聯(lián)合解得

將b＝p1p2/(1+p2),a1＝c2zdz/(1+p2),a2＝c2zc3dz代入上式然后進行中心差分得到

上面三式中，k表示第k次迭代，n表示第n個網格，n+1/2表示第n個半網格；在整個計算區(qū)域上，(56)式和(57)式可以寫成三對角矩陣差分方程，與wlp-fdtd方法相比，這種迭代的因式分解的wlp-fdtd方法將大型稀疏矩陣方程的求解轉變成兩個三對角矩陣方程的求解，于是可以使用追趕法，非常簡單的解得整個計算區(qū)域電磁場分量系數(shù)，最后通過(21)式解得觀測點的電磁場分量，而且這種方法添加了迭代的方案，使得其計算精度比因式分裂的wlp-fdtd要高。

本發(fā)明一種一維高精度迭代的磁化等離子體中的實現(xiàn)方法，流程如圖1所示，具體按照以下步驟實施：

步驟1：輸入模型文件，具體為：

計算區(qū)域大小nz，nz為z軸方向上網格數(shù)；空間步長δz；時間步長δt；真空中的電導率σ，磁導率μ0，介電常數(shù)ε0；等離子體碰撞頻率υ，等離子體頻率ωp，電子回旋頻率ωb；等離子體在計算區(qū)域中的位置；吸收邊界層數(shù)npml與相關參數(shù)κzmax，αzmax，σzmax；κzmax取整數(shù)，κzmax取值范圍為[1，60]；αzmax取值范圍為[01)；σzmax/σopt取值范圍為(0，12]，σopt＝(m+1)/150πδz，m取值范圍為[1，20]，δz取值范圍為λ為源的波長；仿真計算時長tf；迭代次數(shù)k，k≥0且為整數(shù)；加權拉蓋爾多項式的階數(shù)q，q≥0且為整數(shù)；時間尺度因子s，s取值范圍為[10⁹，10¹³]；觀測點；場源參數(shù)。

步驟2：初始化參數(shù)和設置參數(shù)，具體為：

初始化的參數(shù)包括：

將整個計算區(qū)域的電磁場分量系數(shù)整個計算區(qū)域的極化電流密度系數(shù)整個計算區(qū)域的電磁場分量系數(shù)的和整個計算區(qū)域的極化電流密度系數(shù)的和整個計算區(qū)域的輔助變量和其中表示輔助變量表示輔助變量拉蓋爾多項式其中全部初始化為零；磁化等離子體參數(shù)(p1,p2,p3)初始化為p1＝0，p2＝0，pml系數(shù)(c1z，c2z，c3)初始化為c1z＝2/ε0s，c2z＝c1z，c3＝ε0/μ0，其中，ε0是空氣中的介電常數(shù)，μ0是空氣中的磁導率。

設置的參數(shù)具體為：

設置帶有cfs因子的sc-pml吸收邊界的參數(shù)σz,κz,αz；具體為：

σz＝σzmax|z-z0|^m/d^m；

κz＝1+(κzmax-1)|z-z0|^m/d^m；

αz＝αzmax；

式中z0為pml層與非pml截面位置，d是pml吸收邊界的厚度，κzmax取整數(shù)，κzmax取值范圍為[1，60]；αzmax取值范圍為[0，1)；σzmax根據σopt來設置，σzmax/σopt取值范圍為(0，12]；σopt＝(m+1)/150πδz，m取值范圍為[1，20]，δz取值范圍為λ為源的波長；設置pml系數(shù)c1z，c2z；具體為：

c1z＝1/(κzαz+σz+0.5κzε0s)，c2z＝(2αz/ε0s+1)c1z；

設置等離子體參數(shù)為：

步驟3：添加場源到x方向上的電場分量系數(shù)中，設置磁場分量系數(shù)和分別記為初始場值和

其中，所添加的場源的表達式為：

jx(t)＝exp(-(t-t0)²/τ²)

其中，t表示時間，單位為秒；t0,τ為場源參數(shù)。

步驟4：更新計算整個計算區(qū)域的x方向上電場分量系數(shù)具體為：

步驟4.1：電場分量系數(shù)在計算區(qū)域的方程為：

其中，k表示迭代次數(shù)，n表示z軸上第n個計算網格處的位置，n+1/2表示z軸上第n個半網格處的位置，c2z|n表示系數(shù)c2z在z軸方向上第n個網格處的值，c2z|n-1/2表示系數(shù)c2z在z軸方向上第n-1個半網格處的值，δz|n表示z軸方向上第n個網格的大小，第表示第k+1次迭代第n-1個網格處的電場分量系數(shù)的值，

步驟4.2：使用追趕法求解步驟4.1的方程，得到整個計算區(qū)域的電場分量系數(shù)

步驟5：更新計算整個計算區(qū)域的y方向上電場分量系數(shù)具體為：

步驟5.1：電場分量系數(shù)在計算區(qū)域的方程為：

其中，表示第k+1次迭代第n-1個網格處的電場分量系數(shù)的值，

步驟5.2：使用追趕法求解步驟5.1的方程，得到整個計算區(qū)域的y方向上電場分量系數(shù)

步驟6：更新計算整個計算區(qū)域的磁場分量系數(shù)更新方程具體為：

其中

步驟7：更新計算整個計算區(qū)域的磁場分量系數(shù)更新方程具體為：

其中

步驟8：將k+1賦值給k，并判斷迭代次數(shù)k是否達到預設值，若未達到預設值，則返回步驟4，若達到預設值，則執(zhí)行步驟9；

步驟9：更新計算整個計算區(qū)域的極化電流密度系數(shù)更新公式具體為：

步驟10：更新計算整個計算區(qū)域的電磁場分量系數(shù)的輔助變量，更新公式具體為：

其中fζz代表exz,eyz,hxz,hyz，因此表示輔助變量表示電磁場分量系數(shù)的和對z的求導，表示輔助變量的和，表示輔助變量

步驟11：更新計算觀測點處的電磁場分量，更新公式具體為：

上式中u表示電磁場分量ex,ey,hx,hy，r表示電磁場分量的位置，u^q表示q階電磁場分量系數(shù)，是q階加權拉蓋爾多項式，是帶有時間尺度因子s>0的擴展時間，是q階拉蓋爾多項式。

步驟12：將q+1賦值給q，并判斷拉蓋爾多項式的階數(shù)q是否達到預設值，若未達到預設值，則返回步驟3，若達到預設值，則結束。

下面通過實驗對本發(fā)明的效果進行說明：

為了檢驗本發(fā)明方法的正確性、高效性以及高精度，我們計算了9mm厚碰撞磁化等離子體對垂直入射電磁波的反射和透射系數(shù)。入射電磁波為高斯脈沖，加到ex上，其表達式為exp(-(t-t0)²/τ²)，式中t0＝20ps，τ＝5ps。計算區(qū)域為350個網格，每個網格大小為75微米，磁化等離子體占據第201到320的網格，其參數(shù)為ωp＝(2π)·50×10⁹rad/s，ωb＝3.0×10¹¹rad/s，υ＝2.0×10¹⁰hz，其余為真空。完全匹配層吸收邊界放置在計算區(qū)域的兩端，均為10個網格。時間尺度因子為s＝1.885×10¹²，時間步長為0.25ps，仿真時間為tf＝1ns，階數(shù)為200，迭代次數(shù)k為2。pml吸收邊界參數(shù)κzmax＝1，σzmax＝σopt，αzmax＝0。采用本發(fā)明方法、wlp-fdfd方法和因式分裂的wlp-fdtd方法以及解析解計算反射系數(shù)振幅和反射系數(shù)相位，計算結果參見圖2和圖3。從圖中可見，本發(fā)明方法和wlp-fdtd方法、解析解計算結果一致，驗證了本發(fā)明方法的正確性，且本發(fā)明方法的計算精度相比于因式分裂的wlp-fdtd方法的計算精度要高。在運行計算時間上，wlp-fdtd方法需要0.34秒，因式分裂的wlp-fdtd方法需要0.15秒，本發(fā)明專利中的計算方法需要0.23秒，由此可知本發(fā)明專利中的計算方法相對于wlp-fdtd方法計算速度得到了提高，雖然計算速度比因式分裂的wlp-fdtd方法要稍慢點但精度比因式分裂的wlp-fdtd方法要高，尤其是在高頻段，精度提高更明顯。