楊利霞 劉超 李清亮 閆玉波

1) (安徽大學(xué)電子信息工程學(xué)院,合肥 230031)

2) (安徽大學(xué),信息材料與智能感知安徽省實(shí)驗(yàn)室,合肥 230031)

3) (中國(guó)電波傳播研究所,電波環(huán)境特性及模化技術(shù)重點(diǎn)實(shí)驗(yàn)室,青島 266107)

1 引言

電離層作為一種時(shí)變的非均勻磁化等離子體,會(huì)與入射到其中的電磁波發(fā)生復(fù)雜的相互作用[1].一般地,利用高功率尋常波(O 波)對(duì)電離層進(jìn)行加熱實(shí)驗(yàn),由于參量衰變不穩(wěn)定性,泵波可在反射高度附近激發(fā)出電子Langmuir 波和離子聲波[2],同時(shí)波粒相互作用導(dǎo)致O 波轉(zhuǎn)換為Z 波并向電離層更高區(qū)域傳播,進(jìn)一步形成Langmuir 擾動(dòng)[3],Langmuir 擾動(dòng)的形成可引起非相干散射雷達(dá)背向散射功率增強(qiáng),即離子線與等離子體線的增強(qiáng)現(xiàn)象.

電離層加熱是近年來電波傳播領(lǐng)域的熱點(diǎn)問題,電離層局部加熱性導(dǎo)致等離子體的非線性變化[4],可以極大地改變局部地區(qū)的電磁波傳播特性,具有重要的研究?jī)r(jià)值和軍事研究意義.而Langmuir擾動(dòng)是電離層加熱效應(yīng)的重要研究方向之一[5,6],國(guó)外研究人員做了大量的Zakharov 模擬實(shí)驗(yàn)來研究電離層中強(qiáng)擾動(dòng)和弱擾動(dòng)的特性,并將得到的能量譜與觀察到的受激電磁輻射以及雷達(dá)觀測(cè)結(jié)果進(jìn)行比較[7,8].Zakharov 模型對(duì)時(shí)間尺度進(jìn)行了劃分,得到高頻電磁場(chǎng)的復(fù)數(shù)包絡(luò)和等離子體低頻響應(yīng)的方程[9].Mj?lhus 等[10]建立了一個(gè)兩級(jí)模型來表示Langmuir 擾動(dòng)和電磁波之間的耦合,對(duì)Langmuir 波使用Zakharov 模型,對(duì)電磁波使用格林函數(shù)法.Eliasson 等[11,12]進(jìn)行了一系列的研究,對(duì)電離層Langmuir 擾動(dòng)建立了一維(1D)廣義Zakharov 數(shù)值模型.目前國(guó)內(nèi)對(duì)Langmuir 擾動(dòng)本身的研究尚處于起步階段,如武漢大學(xué)的劉默然等[13]、何昉等[14]以及江蘇大學(xué)的朱婷[15]都對(duì)電離層垂直加熱下Langmuir 的擾動(dòng)進(jìn)行了數(shù)值模擬研究.

對(duì)于斜入射[16]電離層加熱下Langmuir 的擾動(dòng)目前國(guó)內(nèi)研究較少,為此,本文基于廣義Zakharov模型,通過將二維(2D)麥克斯韋方程等價(jià)地轉(zhuǎn)換為一維麥克斯韋方程[17,18],避免了用傳統(tǒng)的二維時(shí)域有限差分方法分析該散射問題,這樣將節(jié)省大量?jī)?nèi)存,極大地提高計(jì)算效率,有效降低編程的復(fù)雜程度.同時(shí),深入細(xì)致地研究計(jì)算背景為等離子體的斜入射問題,從而進(jìn)一步研究基于電離層非線性等離子體擾動(dòng)的電磁波傳播特性的全尺度時(shí)域計(jì)算方法,并基于上述方法分析電離層擾動(dòng)的電磁場(chǎng)分布和傳播特性.

2 斜入射非線性等離子體下Langmuir擾動(dòng)分析

2.1 斜入射非磁化等離子體的迭代式推導(dǎo)

TEz斜入射非磁化等離子體的示意圖如圖1所示,其中介質(zhì)區(qū)間1 和3 表示真空,介質(zhì)區(qū)間2表示電離層等離子體.這里選取電磁波的波數(shù)k與z軸夾角為θ,T Ez波在區(qū)間1 中以θ角度斜入射到等離子體區(qū)域中.

圖1 TEz 波斜入射電離層等離子體的傳播模型Fig.1.Propagation model of TEz wave obliquely incident ionospheric plasma.

由圖1 可知,在斜入射情況下,電磁波在電離層中的傳播過程是一個(gè)二維問題.從宏觀角度來看,將麥克斯韋旋度方程在直角坐標(biāo)系下展開,可得 T Ez波在區(qū)間1 的時(shí)域麥克斯韋方程:

根據(jù)文獻(xiàn)[19]中的推導(dǎo)方法,最終(1)式可化為時(shí)域里等價(jià)的一維形式:

式中k1為區(qū)間1 (真空)中的電磁波波數(shù),k1z為k1在z軸方向上的投影,即k1z=k1cosθ.

同理,在區(qū)間2 (等離子體)中,一維斜入射情況下電磁波同樣滿足(2)式形式的一維解,并轉(zhuǎn)換到頻域:

式中k2為區(qū)間2 (電離層等離子體)中的電磁波波數(shù),k2z為k2在z軸方向上的投影,即k2z=k2cos?.

為了解決區(qū)間1,3 與區(qū)間2 的電磁波傳播過程中遇到的相位及邊界問題,由邊界條件和相位匹配原理知將其代入(3)式得

對(duì)方程(4)中的第2 式進(jìn)行整理,可得

式中磁化等離子體的相對(duì)介電常數(shù)為矩陣形式,可寫為

其中σ表示電導(dǎo)率,將(6)式代入(5)式并通過逆傅里葉變換到時(shí)域,可得

在等離子體介質(zhì)中電流密度為

其中,下標(biāo)“ i ”表示離子,下標(biāo)“e”表示電子.

由于離子的質(zhì)量遠(yuǎn)大于電子的質(zhì)量,離子成份的流體速度遠(yuǎn)小于電子成份,離子可以被視為是靜止的.忽略離子部分后將(8)式代入(7)式,最后聯(lián)立(4)式的第一個(gè)式子并進(jìn)行逆傅里葉變換,可得

同理推導(dǎo)Ey1D,Hx1D具有相似的形式:

電離層等離子體可被看成是由電子流體和離子流體組成的導(dǎo)電的連續(xù)介質(zhì),因此描述等離子體運(yùn)動(dòng)采用流體力學(xué)方程.在雙流體力學(xué)描述中,帶電粒子的數(shù)密度和流體速度分別滿足連續(xù)性方程和動(dòng)量方程:

式中nα為粒子數(shù)密度,υα為粒子的流體速度,qα表示粒子電荷;磁感應(yīng)強(qiáng)度B=B0+B1,其中B0是地磁場(chǎng),B1是波磁場(chǎng);mα是粒子質(zhì)量;Pα是壓強(qiáng);γα是與粒子間碰撞等因素有關(guān)的阻尼算子.

將(11)式中的電子數(shù)密度連續(xù)性方程,經(jīng)過高斯定理推導(dǎo)出如下形式:

對(duì)于高頻波,(11)式中的電子動(dòng)量方程變?yōu)?/p>

(9)式、(10)式、(12)式和(13)式構(gòu)成了高頻電磁波斜入射電離層等離子體的非線性耦合方程組,也是此研究的理論基礎(chǔ).

時(shí)域有限差分方法可以直接求出Maxwell 方程組顯式下的離散形式,(9)式和(10)式進(jìn)行離散后得到磁場(chǎng)強(qiáng)度的迭代方程為

從迭代式(14)—(16)可以看出,電場(chǎng)的垂直分量E⊥與磁場(chǎng)的垂直分量以及等離子體密度、流體速度密切相關(guān),而電場(chǎng)的平行分量Ez僅與其本身、電子的數(shù)密度和z方向的流體速度有關(guān).

對(duì)(12)式關(guān)于空間的一階偏導(dǎo)數(shù)選取后向差分近似,得到電子密度的離散迭代方程為

對(duì)(13)式等號(hào)右側(cè)的第一項(xiàng)作近似處理,再通過矩陣求解得到電子流體速度的離散迭代方程為

電磁波在電離層中斜入射傳播過程的時(shí)域數(shù)值完整迭代步驟如圖2 所示.

圖2 FDTD 電磁場(chǎng)時(shí)域的逐步迭代流程圖Fig.2.Step by step iterative flow chart of FDTD electromagnetic field in time domain.

2.2 數(shù)值算法的穩(wěn)定性條件

1)時(shí)間離散間隔的穩(wěn)定性要求

在使用FDTD 算法對(duì)電磁波在電離層中的傳播進(jìn)行數(shù)值仿真時(shí),需要考慮等離子體與電磁波相互作用的影響.等離子體中空間和時(shí)間離散間隔之間應(yīng)當(dāng)滿足的充分條件[20]為

其中wpe為等離子體特征頻率.(19)式表明一維情況下,時(shí)間間隔必須等于或小于波以光速通過Yee元胞所需的時(shí)間.

2)數(shù)值色散對(duì)空間離散間隔的穩(wěn)定性要求

有限差分代替波動(dòng)方程的二階導(dǎo)數(shù)過程中,離散處理將會(huì)導(dǎo)致波的色散,此時(shí)需要對(duì)空間離散間隔做出要求來降低色散,一維情況下要求

式中λ對(duì)應(yīng)介質(zhì)中的最小波長(zhǎng).

本次數(shù)值模擬選取λ=60 m ,Δz=1 m,并使Δt滿足Courant 條件,本文在50—400 km 全范圍內(nèi)以1 m 步進(jìn)計(jì)算.

2.3 數(shù)值模型

本文建立的電磁波斜入射電離層的數(shù)學(xué)模型如圖3.數(shù)值模型采取的實(shí)際參數(shù)來自EISCAT 在挪威的電離層加熱實(shí)驗(yàn)[21,22],考慮到參量不穩(wěn)定性發(fā)生的反射區(qū)域高度范圍一般在200—340 km,假設(shè)沿z軸垂直向上為正方向,水平方向正北為y軸正方向,水平向西為x軸正方向.計(jì)算區(qū)域的范圍在垂直方向延伸400 km.為了降低計(jì)算量對(duì)計(jì)算機(jī)內(nèi)存的要求,本文的數(shù)值模擬將電離層剖面的尺寸壓縮為實(shí)際尺寸的1/20.

圖3 電磁波加熱電離層傳播模型Fig.3.Ionospheric propagation model heated by electromagnetic wave.

2.4 數(shù)值仿真參量

在z=50 km 處通過連接邊界,采用線極化的調(diào)制高斯脈沖作為入射源,Ex=1.5 sin(ω0t)exp[-(z-50×103)/108]V/m.電磁波在電離層的傳播過程中發(fā)生模式轉(zhuǎn)換進(jìn)而產(chǎn)生Langmuir 擾動(dòng)現(xiàn)象,必須滿足兩個(gè)條件:1)入射波頻率應(yīng)當(dāng)略小于 F2層的最大臨界頻率[23,24];2)在臨界層內(nèi)存在等離子體密度或者電場(chǎng)強(qiáng)度的擾動(dòng)[25].根據(jù)所給參數(shù),計(jì)算出 F2層最大截止頻率fmax≈6.33 MHz,這里選取加熱電磁波頻率f0=5 MHz 以滿足第1 個(gè)條件.第2 個(gè)條件是電磁波對(duì)電離層加熱引起了參量的不穩(wěn)定性,會(huì)使等離子體中存在靜電波擾動(dòng).

在模擬計(jì)算區(qū)域中,上下邊界的截?cái)嗵幉捎靡痪S的Mur 吸收邊界.根據(jù)吸收邊界原理和連接邊界理論[26],在斜入射情況下空間步長(zhǎng)改為 Δz在z方向上的投影,具體修正為如下形式:

左截?cái)噙吔缣?/p>

TEz波只需在原有的連接邊界形式上乘以角度因子 c osθ,具體修正為如下形式:

2.5 O 波、X 波和Z 波

對(duì)于高頻電磁波,考慮垂直傳播的波,假設(shè)電磁波傳播方向k沿E1×B1方向,因此電場(chǎng)可能有兩種基本方向(如圖4 所示):E1//B0和E1⊥B0.當(dāng)E1平行于磁場(chǎng)B0時(shí),這種波被稱為尋常(O)波[27],得到O 波的色散關(guān)系為結(jié)果表明尋常波的傳播不受磁場(chǎng)影響.對(duì)于E1垂直于磁場(chǎng)B0的情況,電子響應(yīng)隨電場(chǎng)而運(yùn)動(dòng),但洛倫茲力的作用使電子不能沿固定方向運(yùn)動(dòng),而要產(chǎn)生另一方向的速度分量.為了區(qū)別于和B0無關(guān)的尋常波,把和B0有關(guān)的波動(dòng)稱為非尋常(X)波,非尋常波是由部分橫波和部分縱波構(gòu)成的電磁波.非尋常波的色散方程為

圖4 O 波與X 波Fig.4.O wave and X wave.

其中n為折射率,等離子體對(duì)于擾動(dòng)的響應(yīng)是一個(gè)頻率為的振蕩,該頻率稱為高混雜頻率.

Z 波模式是特殊(X)模式的低頻分支的空間物理符號(hào),它是等離子體的一個(gè)內(nèi)?；蚍@模,其頻率限定在右旋截止頻率WL和高混雜頻率WUH之間[28].

3 模擬結(jié)果

3.1 截止和共振

為了分析波的傳播特性,引入兩個(gè)重要的概念:截止和共振[29].當(dāng)波的折射率變?yōu)榱慊虿〝?shù)趨于零時(shí),在等離子體中出現(xiàn)波的截止.由O 波的色散關(guān)系式知,波數(shù)趨于0 時(shí)其截止頻率等于電子等離子體頻率wpe,當(dāng)w＞wpe時(shí),O 波可以在等離子體中傳播.由(24)式得,X 波在n2=0 ,k=0 情況下截止,因?yàn)閚2＜0 使得波傳播因子變?yōu)檎穹p因子,意味著波在介質(zhì)中傳播時(shí)很快衰減,最終被截 止.而在n2→∞,k→∞條件下w與k無關(guān),這樣相速度、群速度都為0,波不能傳播,出現(xiàn)共振.在截止點(diǎn)和共振點(diǎn),波都不能傳播,一般來說波在截止點(diǎn)被反射,在共振點(diǎn)被吸收.

令n2=0 求解(24)式,w方程應(yīng)該有4 個(gè)根,其中只有兩個(gè)根是合理的(w＞ 0),它們的截止頻率分別稱為左旋、右旋截止頻率,如圖5 所示.

圖5 w-k 色散曲線Fig.5.w-k dispersion curve.

通過對(duì)O 波和X 波的色散分析可知,對(duì)于O 波,當(dāng) 0＜w ＜wpe時(shí)不能傳播,w＞wpe時(shí)才能傳播.對(duì)于X 波,在 0＜w ＜wL和wUH＜w ＜wR兩 個(gè)頻率區(qū)域不能傳播,不能傳播的頻區(qū)稱為截止帶;當(dāng)wL＜w ＜wUH和w＞wR時(shí),X 波可以傳播,對(duì)應(yīng)的頻區(qū)稱為傳播帶.

3.2 斜入射下的電磁波傳播特性

通過數(shù)值仿真模擬斜入射角為7°時(shí)的電離層傳播情況,結(jié)果如圖6 所示.本文數(shù)值模型的計(jì)算是在激勵(lì)源(50 km,圖6 粉色圓點(diǎn)處)加熱一段時(shí)間后開始進(jìn)行,電場(chǎng)的初始值在下邊界是一個(gè)沿著x方向線性極化的波,在進(jìn)入電離層等離子體區(qū)域前先經(jīng)過一段中性大氣(50—85 km,紅色虛線區(qū)域),該傳播過程中包絡(luò)形狀不變.

圖6 背景電子密度及t=0.249 ms 時(shí)刻電場(chǎng) E x 分量幅值分布Fig.6.Background electron density and amplitude distribution of electric field E x component at t=0.249 ms.

入射波在電離層等離子體中傳播一段距離后出現(xiàn)了y方向上的分量,如圖7 所示,這是由法拉第旋轉(zhuǎn)現(xiàn)象引起的,并且電場(chǎng)y分量的包絡(luò)形狀與入射波基本一致.

圖7 背景電 子密度 及t=0.821 ms 時(shí)刻電 場(chǎng) E x 和 E y 分量幅值分布Fig.7.Background electron density and distribution of electric field E x and E y component amplitude at t=0.821 ms.

如圖8 所示,隨著等離子體密度的增加,電磁波的波形開始發(fā)生變化,X 波反射點(diǎn)大約在271 km處(紅點(diǎn)所示),電磁波逐步接近X 波反射點(diǎn)的過程中,會(huì)導(dǎo)致波包出現(xiàn)分界點(diǎn).

圖8 背景電 子密度 及t=0.90 ms 時(shí)刻電 場(chǎng) E x 和 E y 分量幅值分布Fig.8.Background electron density and distribution of electric field E x and E y component amplitude at t=0.90 ms.

在達(dá)到X 波反射點(diǎn)后還有O 波繼續(xù)向前傳播,O 波在接近其反射點(diǎn)的過程中逐漸從左旋圓極化變成沿地磁場(chǎng)方向的線極化,圖9 所示的O 波反射點(diǎn)附近激發(fā)出了明顯的z向局部振蕩的電場(chǎng)分量.從圖9 可以明顯看出,大約在278 km 附近有明顯的截止特征,該處為O 波的反射點(diǎn),與根據(jù)Appleton-Hartree 公式計(jì)算所得的O 波反射高度相符.

圖9 背景電 子密度 及t=1.03 ms 時(shí)刻電 場(chǎng) E x ,E y 和Ez分量幅值分布Fig.9.Background electron density and distribution of electric field E x ,E y and E z component amplitude at t=1.03 ms.

由圖10 可以看出,z向分量電場(chǎng)大約在1.2 ms之后到達(dá)O 波反射區(qū)域,并形成駐波.圖11 是在O 波轉(zhuǎn)換區(qū)域(277—280 km)以及更高區(qū)域(310—320 km)隨時(shí)間變化的靜電擾動(dòng)及電子密度(ns=ne-ni0)圖像.

圖10 背景電子密度及t=1.27 ms 時(shí)刻電場(chǎng) E x ,E y 和Ez分量幅值分布Fig.10.Background electron density and distribution of electric field E x ,E y and E z component amplitude at t=1.27 ms.

從圖11(a)可以看出,大約在1.2 ms,有參量不穩(wěn)定性引起的靜電波擾動(dòng)導(dǎo)致O 波轉(zhuǎn)換為Z 波,并向電離層的更高區(qū)域傳播.圖9 中電場(chǎng)在O 波反射點(diǎn)附近的快速振蕩也是由于參量不穩(wěn)定性過程中的Langmuir 波造成的,即高頻Langmuir波的連續(xù)振蕩引起電子運(yùn)動(dòng).如圖11(b)所示,在相互作用區(qū)域周圍越來越多的電子參與波的相互作用.從圖11(c)可以明顯看到O 波以上,位于314—316 km 之間依然有波的存在,此處為Z 波的轉(zhuǎn)換區(qū)域[30].從圖11(c)可以看出在Z 波反射區(qū)域,其電場(chǎng)分量的高度非常窄,這里是Z 波經(jīng)過參量不穩(wěn)定性激發(fā)出的短波段Langmuir 波,對(duì)應(yīng)的電子密度如圖11(d)所示.

圖11 O 波轉(zhuǎn)換 區(qū)域的 靜電擾動(dòng)及 ns 擾 動(dòng) (a) O 波轉(zhuǎn)換區(qū)域 | Ez| 變 化;(b) O 波轉(zhuǎn)換區(qū)域 | ns| 變 化;(c) Z 波轉(zhuǎn)換區(qū)域 | Ez| 變 化;(d) Z 波轉(zhuǎn)換區(qū)域 | ns| 變 化Fig.11.Electrostatic disturbance and ns disturbance in O-wave conversion region:(a) Variation of O-wave conversion region with respect to | Ez| ; (b) variation of O-wave conversion region with respect to | ns| ;(c) variation of Z-wave conversion region with respect to | Ez| ; (d) variation of Z-wave conversion region with respect to | ns| .

4 結(jié)論

通過二維麥克斯韋方程等價(jià)轉(zhuǎn)換成一維麥克斯韋方程組,與雙流體力學(xué)方程進(jìn)行聯(lián)立,建立了斜入射電離層等離子體的數(shù)值模型,模擬了大功率電磁波在較小傾斜角度加熱電離層情況下O 波的參量不穩(wěn)定性.結(jié)果表明:在毫秒量級(jí)的時(shí)間尺度內(nèi),大功率高頻電磁波在小角度的斜入射加熱下,O 波反射區(qū)域有參量不穩(wěn)定性引起的靜電波擾動(dòng),O 波發(fā)生明顯的衰減并激發(fā)了Langmuir 波的形成以及Langmuir 擾動(dòng)所引起的O 波轉(zhuǎn)換為Z 波現(xiàn)象.其中,Z 波在1.2 ms 左右到達(dá)O 波反射區(qū)域并形成駐波.