黏性各向異性磁流體Kelvin-Helmholtz不穩定性: 二維數值研究*

劉迎 陳志華 鄭純

(南京理工大學, 瞬態物理國家重點實驗室, 南京 210094)

(2018 年 9 月 21 日收到; 2018 年 11 月 19 日收到修改稿)

利用 corner transport upwind 和 constrained transport算法求解非理想磁流體動力學方程組, 對勻強平行磁場作用下, 黏性各向異性等離子體自由剪切層中的Kelvin-Helmholtz不穩定性進行了數值模擬. 從流動結構、渦結構演化、磁場分布、橫向磁壓力、抗彎磁張力等角度對各向同性和各向異性黏性算例結果進行了討論, 分析了黏性各向異性對Kelvin-Helmholtz不穩定性的影響. 結果表明, 黏性各向異性比黏性各向同性更利于流動的穩定. 其穩定性作用是由于磁感線方向上剪切速率降低導致界面卷起程度和圈數的降低, 并使卷起結構中小渦產生增殖、合并, 破壞了渦的常規增長, 從而導致流動的穩定. 黏性各向異性對橫向磁壓力的影響比對抗彎磁張力更大.

1 引 言

速度剪切層中因擾動而產生卷起、合并的Kelvin-Helmholtz (KH)不穩定性廣泛存在于自然界中, 它對于研究行星的演化[1,2]、地磁場變化[3?5]、核聚變控制[6?8]、氣象變化[9, 10]等具有重要意義. 早期的工作大多集中于各向同性流體[11?16]. 然而在等離子體中, 粒子的碰撞和輸運會對整體的宏觀物理屬性產生較大影響. 由于洛倫茲力同時垂直于帶電粒子的速度方向和磁場方向, 因此等離子體在平行于磁場方向和垂直方向上的物理屬性有著明顯的差異. 通常, 平行于磁場方向上的輸運系數遠大于垂直于磁場方向上的輸運系數, 這是實際中存在的各向異性. 考慮到各向異性等離子體在實際中大量存在, 因此對于各向異性等離子體KH不穩定性的研究, 尤其是各向異性對流動結構和演變影響以及對不穩定性影響的研究十分必要.

前人關于各向異性磁流體KH不穩定性的研究不多. 早在 1965 年, Talwar[17]研究了均勻磁場作用下, 正切速度間斷分布的壓力各向異性等離子體流動問題. 結果發現, 對于剪切速度在一定范圍內, 單調不穩定性是可能的. 1975 年, Duhau 和Gratton[18]研究了可壓縮壓力各向異性等離子體在剪切磁場和流向均勻磁場作用下的KH不穩定性問題, 發現壓力各向異性在穩定性方面起著很重要的作用. 隨后, Srivastava 和 Vyas[19]同樣研究了沿直磁場方向剪切流動的非均勻(流速和等離子體密度都沿磁場垂直方向上變化)各向異性等離子體中的KH不穩定性問題, 討論了某些極限情況下回轉黏度與不均勻性對不穩定性的影響.

1985 年, Choudhury 和 Patel[20]使 用 雙 絕 熱Chew-Goldberger-Low (CGL)方程研究了線性速度分布的各向異性等離子體中的KH不穩定性問題, 得到不同馬赫數-波數平面內擾動增長率結果.隨后, Choudhury[21]將研究推進到正切速度分布的各向異性等離子體情形. 1996 年, Ruderman 等[22]研究了具有各向異性黏性及熱導率的可壓縮等離子體在正切流速分布情況下的流動問題, 并討論了不可壓的極限情形. 他們發現, 熱導率在不可壓情形下作用消失, 此時只有黏性有影響; 當剪切層兩端黏性不同時, 黏性有失穩作用; 黏性的存在將會降低無黏情形下的臨界剪切速度. 2002年,Brown和Choudhury[23]使用線性不穩定性分析理論及廣義多方定律 (generalized polytrope laws),研究了具有正切間斷速度分布的壓力各向同性和壓力各向異性磁流體中的KH不穩定性問題, 繪制了逆磁壓比-擾動傳播角度平面上的穩定區域, 闡明了磁場與壓縮性的影響. 2010 年, Prajapati和Chhajlani[24]利用廣義多方定律, 研究了流向磁場作用下等離子體中壓力各向異性和流速KH不穩定性的影響, 使用了雙絕熱CGL方程組以及磁流體力學(MHD)方程組兩種模型討論流向擾動的穩定性. 結果發現, 與 CGL 方程組相比, MHD 方程組的增長率更大; 壓力各向異性具有破壞穩定的作用, 磁場對系統有穩定作用.

前人關于各向異性大多基于黏性分層均勻分布, 或壓力各向異性, 少有研究黏性各向異性問題.本文使用數值模擬方法求解非理想MHD方程組,研究黏性各向異性對KH不穩定性的作用與機理.在第2部分給出研究的物理模型與模擬所用的數值方法, 使用的是 CTU + CT (corner transport upwind + constrained transport)算 法[25], 其 中CTU算法用于多維積分, CT算法用于保證磁場散度為零. 在第3部分中給出數值模擬結果及討論;最后是結論.

2 模型與計算方法

忽略霍爾效應、雙極擴散效應和相對論效應,帶有各向異性黏性項的非理想MHD方程組的守恒形式為

上述方程組中用到了磁場無散度約束條件?·B=0. 其中流體的總能體積密度為

各向同性剪切應力張量為

各向異性剪切應力張量為

上述方程中,ρ為流體密度,u為流體速度矢量,p為等離子體壓力,I為二階單位張量,B是磁感應強度,b為磁場單位矢量,μm為磁導率,μ為動力黏性系數 (下標“0”表示各向同性分量, 下標“||”表示各向異性分量). 狀態方程采用理想氣體狀態方程, 絕熱指數γ=1.4 .

本文使用有限體積法對方程進行離散, 采用PPM (piecewise parabolic method)三階空間重構和 HLLD (HLL-discontinuties)通量格式. 使用無方向分裂的CTU算法求解上述方程組, 應用CT算法處理磁場無散度約束.

為方便研究, 本文選擇二維平面剪切層為研究對象, 如圖1所示. 考慮兩股在x方向上相向流動的自由剪切層流動, 假設流體是理想導電、具有各向異性黏性, 水平方向上流速分布為雙曲正切間斷,

上下層速度差 ?U=0.5a, 這里a為聲速. 整個流場附加大小為B0的流向均勻磁場,

而y方向(縱向)上附加高斯衰減的正弦速度擾動,

其中ε為振幅,σ為波幅度衰減率,λ為波長,kx=2π/λ是x方向上的擾動波數. 為防止擾動對流動干擾過大, 這里選擇小振幅的擾動ε=0.03 , 并且在過渡層內快速衰減σ=2δ.

初始時刻, 來流中上、下層具有相同的初始壓力p0、密度ρ0、比熱比γ=1.4 、聲速以及阿爾文速度為了獲得渦結構最大分辨率, 設置擾動波長λ=L. 為了使得不穩定性增長較為明顯, 選擇相對較窄的過渡層δ=L/50. 本文中, 采用L為單位長度, 聲音跨越長度L的時間L/a為單位時間,作為單位磁場強度以無量綱化計算結果.

并將其與MHD方程組耦合求解.

圖1 計算模型Fig.1. Sketch of the computational model.

3 結果與討論

圖2為不同時刻剪切層示蹤劑濃度場截圖, 給出了各向同性黏性和各向異性黏性磁流體中KH不穩定性的形成與演變過程. 從圖2可知, 各向同性黏性磁流體剪切層在擾動下界面發生卷起,然而由于磁場的作用, 卷起結構被斜向拉長, 但不能繼續卷起形成貓眼渦結構(圖2(a)). 對于各向異性黏性情形(圖2(b)), 剪切層仍然發生了卷起, 被斜向拉長, 但也不能形成典型的貓眼渦結構, 且卷起結構中心處出現了波紋狀界面. 對比圖2(a)和圖2(b)可知, 圖2(b)中卷起的圈數減少. 這些特征暗示各向異性黏性的加入具有致穩效果.

由于黏性影響渦的擴散, 下面從渦結構演化的角度來分析其作用過程. 最近Liu等[26]提出的Rortex方法表現出較好的渦結構展示能力, 本文使用該方法對渦結構進行可視化. 圖3給出了各向同性黏性算例中Rortex場在不同時刻的截圖. 從圖 3 可以看出,t= 5.5時形成了順時針旋轉的“貓眼”渦(圖3(a)). 隨著界面繼續卷起, 渦結構增大,同時螺旋速度降低, 渦內速度方向開始發生改變(圖3(b)), 這是由于黏性降低了剪切作用所致. 當卷起結構在磁場作用下被斜向拉長, 流速不均勻導致內部出現“乙”形渦流(圖3(c)), 并且渦內部有發生分裂擴散的趨勢.t= 8.0時, 渦結構外部上下邊緣向內凹陷, 而渦流中心已凹陷, 渦中心左右兩側各有兩個一大一小的渦流, 均為順時針方向(圖3(d)).到t= 9.0時, 斜向拉長的渦流外側上下緣各形成一狹長逆時針方向渦流, 左右兩側各分布兩指甲狀順時針渦流, 而中心處逆時針渦流被上下兩個狹長順時針小渦流夾住 (圖 3(e)). 在t= 10 時刻, 左右兩側指甲狀渦以及上下兩側狹長渦圍繞卷起結構中心逆時針旋繞, 它們在內部渦流向外擴散過程中被擠壓成月牙形, 而中間的渦流擴大, 左右兩端在旋轉后增殖出兩個小圓形渦(圖3(f)). 上述過程表明了黏性對剪切層渦結構的分裂增殖以及擴散作用, 使流動中出現了逆時針渦流, 從而感應出逆時針的速度場, 以對抗順時針旋轉的卷起.

圖4為各向異性算例中Rortex方法展示的渦結構演化過程. 其演化與圖3類似, 但細節稍有不同, 下面來對比它們的不同點. 在t= 5.5 時, 圖 3(a)中的渦流中心處(藍色區域)分布均勻, 并且比左右兩側(藍色區域)要深, 這表明渦流是均勻的且中心處旋轉速度最大. 但在圖4(a)中, “貓眼”渦流中心為條紋狀的紋理(藍黃相間), 且中心處的顏色不如左右兩側的藍色區域深, 這說明中心處的流動出現了不均勻現象, 且中心處旋轉速度低于兩側.到t= 7.0 時, 圖 3(b)中心處旋轉速度稍微減小,卻與左右兩側大小幾乎相當, 而圖4(b)中心旋轉速度與左右兩側(深藍色區域)相差較多. 這說明黏性各向異性加入后減弱了卷起結構中心處的渦旋速度. 隨后在t= 7.5 時, 圖 3(c)中流線開始扭曲, 中心處形成順時針旋轉的傾斜“乙”字形渦流,而圖4(c)中心卻分離出一對順時針旋轉的橢圓小渦. 到t= 8.0 時, 圖 3(d)中心兩側出現一大一小渦流, 且中心處流線向內凹陷, 然而圖4(d)中心處兩側只有兩個渦流, 且中心處流線并未明顯向內凹陷. 到t= 9.0 時, 圖 3(d)中間形成逆時針旋轉渦流, 整體結構上下兩側為逆時針旋轉的渦流, 內部還分散著位置對稱的小渦. 但圖4(d)中則結構簡單, 中心處新產生逆時針旋轉渦流, 上下兩側流線發生彎曲, 但尚未閉合形成渦. 這一系列現象表現出各向異性黏性使流動結構變得更加規則, 渦的生成和合并減慢. 到t= 10 時, 圖 3(e)中卷起結構內部中心處, 幾對渦合成一個較大的逆時針旋轉橢圓型渦流, 邊緣處的渦被擠壓成月牙形. 圖4(e)中心也形成了逆時針旋轉的橢圓型長渦, 但上下兩個渦流并未被強烈擠壓.

圖2 不同時刻自由剪切層示蹤劑濃度場 (a)各向同性黏性 (Re = 105); (b) 各向異性黏性 (Re0 = 105, Re|| = 100)Fig.2. Tracer concentration field at different times: (a) Isotropic viscous case (Re = 105); (b) anisotropic viscous case (Re0 = 105,Re|| = 100).

圖3 不同時刻, 各向同性黏性 (Re = 105)算例渦結構云圖及流線 (a) t = 5.5; (b) t = 7.0; (c) t = 7.5; (d) t = 8.0; (e) t = 9.0;(f) t = 10Fig.3. Rortex field with streamlines at different times in isotropic viscous fluid (Re = 105): (a) t = 5.5; (b) t = 7.0; (c) t = 7.5;(d) t = 8.0; (e) t = 9.0; (f) t = 10.

從上述對比中可以發現, 各向異性加入后, 渦流中的渦旋速度降低, 渦的生成和增殖速率降低.由此可見, 流動的穩定程度提高了.

圖4 不同時刻各向異性黏性 (Re0 = 105, Re|| = 100) 算例渦結構云圖及流線 (a) t = 5.5; (b) t = 7.0; (c) t = 7.5; (d) t = 8.0;(e) t = 9.0; (f) t = 10Fig.4. Rortex field with streamlines at different times in anisotropic viscous fluid (Re0 =105, Re|| = 100): (a) t = 5.5; (b) t = 7.0;(c) t = 7.5; (d) t = 8.0; (e) t = 9.0; (f) t = 10.

圖5 不同時刻剪切層中磁場云圖和磁感線 (a) 各向同性黏性 (Re = 105); (b) 各向異性黏性 (Re0 = 105, Re|| = 100)Fig.5. Magnetic field with field lines at different times: (a) Isotropic viscous case (Re = 105); (b) anisotropic viscous case (Re0 =105, Re|| = 100)

圖5為不同時刻各向同性與各向異性黏性算例中磁場強度和磁感線的分布. 先考察各向同性黏性情形 (圖 5(a)). 在t= 5.5 時, 剪切界面發生卷起, 由于磁凍結效應, 磁感線也隨之卷起, 形成螺旋狀結構. 卷起使磁感線被拉伸和彎曲, 導致磁場被放大, 因此圖中能清晰地看到卷起的磁感線上磁場強度約為周圍磁場的 3 倍以上.t= 8.0 時, 界面進一步卷起并被斜向拉伸, 整體結構約45°傾斜,圖中可見磁感線卷起程度提高且在卷起外緣(紅色區)磁場放大程度較強. 然而到t= 10 時, 整個磁感線結構進一步被斜向拉長, 而且有向內被擠壓的趨勢, 因而磁場放大程度較強的(紅色)區域分布較廣. 但在卷起結構內部, 磁感線卷曲程度變得比t= 9.0時要低, 這反映了抗彎磁張力的穩定性作用.

對于各向異性黏性情形(圖5(b)), 可與各向同性黏性情形分別對比. 首先, 在結構方面, 各向異性黏性情形中磁感線卷起結構中的螺旋的圈數少于各向同性黏性情形, 中心處幾乎沒有卷曲形狀. 其次, 在磁場強度方面, 各向異性黏性情形中卷起界面上磁場放大程度低于各向同性黏性情形,且磁場被放大的區域分布范圍不集中. 這表明, 各向異性黏性的加入, 使得沿磁感線方向上的剪切程度降低, 從而降低了磁感線卷曲程度, 導致磁場放大程度的下降, 使得流動變得穩定.

圖6為流場中縱向總動能隨時間的變化. 可見各向同性黏性算例和各向異性黏性算例都經歷了線性增長與飽和后的非線性階段. 然而不同的是,各向異性黏性算例的線性增長率以及飽和水平都比各向同性黏性算例的低, 且線性增長時間稍微增加. 這從動能的角度表明, 各向異性黏性相比于各向同性黏性, 具有促進穩定的作用.

圖6 縱向總動能隨時間的變化Fig.6. Evolution of the longitudinal total kinetic energy Eky .

最近, Liu 等[27]研究發現, 磁場對 KH 不穩定性的致穩作用源于磁能局部放大所產生的橫向磁壓力(fmp)和抗彎磁張力(fmt), 下面分析這些磁場量的變化.

圖7為平均磁能密度相對于初始時刻放大倍數隨時間的變化. 可見, 在t= 3.5 之前平均磁能密度幾乎保持不變, 之后開始非線性上升, 直到大約在t= 9.5時第一次到達峰值. 其中各向同性黏性算例中平均磁能密度在峰值時放大到9倍左右,而各向異性黏性算例中平均磁能密度則放大到6倍左右. 這正好解釋了圖4中局部磁場放大的區別.

圖7 平均磁能密度放大倍數隨時間的變化Fig.7. Evolution of the amplification factor of average magnetic energy.

圖8為平均橫向磁壓力和抗彎磁張力隨時間的變化. 圖8(a)中各向同性情形橫向磁壓力(紅色實線)初期幾乎線性上升到t= 4左右, 此后開始非線性增長, 直到t= 10 左右飽和. 對于各向異性黏性算例(藍色點線), 其曲線變化趨勢與各向同性黏性算例相似, 但增長速度和幅值始終低于各向同性黏性情形. 這是因為各向異性黏性降低界面卷起程度, 減小磁能密度的放大程度, 從而降低了橫向磁壓力. 圖8(b)中的抗彎磁張力則有所不同, 各向同性黏性與各向異性黏性算例對應曲線幾乎重合,二者都在t= 0—0.5左右快速線性上升, 然后非線性緩慢上升, 直到t= 10. 這表明各向異性黏性幾乎對抗彎磁張力無影響.

圖9為不同算例中流場總渦度擬能隨時間的變化. 可見, 曲線在t= 0 到t= 3.5 左右幾乎重合, 且都是緩慢下降. 隨后, 各向同性黏性算例對應的曲線開始上升, 在t= 7之后幾乎直線上升;而各向異性黏性算例仍然沿之前趨勢緩慢下降到t= 7.5 左右才開始緩慢上升. 到t= 10 時, 各向異性黏性算例中渦度擬能只有各向同性的25%左右. 這表明, 黏性和磁場在KH不穩定波的線性增長階段發揮了穩定性作用, 緩慢降低了渦度擬能.隨后在非線性發展階段, 由于黏性的擴散作用, 導致流動中生成了許多細小的渦, 而之前擾動動能則一部分轉化為磁能, 另一部分則轉化為渦度擬能.各向異性黏性則由于其加強了磁感線方向的黏性,使得沿著磁感線方向上的剪切速度變慢, 從而有效降低了渦量, 因此曲線上升速度不及各向同性黏性算例.

圖8 平均橫向磁壓力和抗彎磁張力隨時間的變化 (a) 橫向磁壓力; (b) 抗彎磁張力Fig.8. Evolution of the average transverse magnetic pressure and anti-bending magnetic tension: (a) Transverse magnetic pressure; (b) anti-bending magnetic tension.

圖9 流場總渦度擬能隨時間的變化Fig.9. Evolution of the total enstrophy.

4 結 論

利用CTU + CT算法對非理想磁流體方程組進行求解, 模擬了勻強平行磁場作用下, 具有各向異性黏性剪切層中的KH不穩定性, 研究了黏性各向異性對KH不穩定性的影響和作用機理.

模擬結果表明, 各向異性黏性具有穩定作用,在同樣情況下, 其致穩效果比各向同性黏性稍強.其穩定作用是因為: 各向異性黏性中沿磁感線方向上的黏性系數比垂直于磁感線方向上的黏性系數要大得多, 導致界面卷起過程中沿著剪切面上的剪切速度變慢, 從而使卷起結構內部的卷曲圈數降低, 并且卷曲程度下降, 最終流動不再繼續往不穩定方向發展.

在該過程中, 黏性使得卷起結構中小渦產生增殖, 特別是在邊緣和內部分裂出數個逆時針的小渦. 這些小渦還會發生合并, 破壞了卷起結構常規增長, 從而使得流動變得穩定. 各向異性黏性對磁場的影響主要體現在降低橫向磁壓力方面, 對抗彎磁張力的影響較小.