摘要:利用corner transport upwind和constrained transport算法求解非理想磁流体动力学方程组, 对匀强平行磁场作用下, 黏性各向异性等离子体自由剪切层中的Kelvin-Helmholtz不稳定性进行了数值模拟. 从流动结构、涡结构演化、磁场分布、横向磁压力、抗弯磁张力等角度对各向同性和各向异性黏性算例结果进行了讨论, 分析了黏性各向异性对Kelvin-Helmholtz不稳定性的影响. 结果表明, 黏性各向异性比黏性各向同性更利于流动的稳定. 其稳定性作用是由于磁感线方向上剪切速率降低导致界面卷起程度和圈数的降低, 并使卷起结构中小涡产生增殖、合并, 破坏了涡的常规增长, 从而导致流动的稳定. 黏性各向异性对横向磁压力的影响比对抗弯磁张力更大.
关键词: 磁流体动力学/
开尔文-亥姆霍兹不稳定性/
黏性各向异性/
涡

English Abstract

全文HTML

--> --> -->

速度剪切层中因扰动而产生卷起、合并的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部分中给出数值模拟结果及讨论; 最后是结论.

忽略霍尔效应、双极扩散效应和相对论效应, 带有各向异性黏性项的非理想MHD方程组的守恒形式为

$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho {{u}}} \right) = 0 , $

$ \frac{{\partial \rho {{u}}}}{{\partial t}} + \nabla \cdot \left( {\rho {{uu}} + p{{I}} + \frac{{{{B}} \cdot {{B}}}}{{2{\mu _{\rm m}}}} - \frac{{{{BB}}}}{{{\mu _{\rm m}}}}} \right) = - \nabla \cdot {{\tau }} , $

$ \begin{split}& \frac{{\partial \rho {e_{\rm{t}}}}}{{\partial t}} + \nabla \cdot \left( {\rho {e_{\rm t}}{{u}} + p{{u}} + \frac{{{{B}} \cdot {{B}}}}{{2{\mu _{\rm m}}}}{{u}} - \frac{{{{BB}} \cdot {{u}}}}{{{\mu _{\rm m}}}}} \right)\!\!\!\!\!\! \\=& - \nabla \cdot ({{\tau }} \cdot {{u}}) , \end{split}$

$ \frac{{\partial {{B}}}}{{\partial t}} = \nabla \times ({{u}} \times {{B}}) .$

上述方程组中用到了磁场无散度约束条件$\nabla \cdot {{B}}=0$. 其中流体的总能体积密度为

$ \rho {e_{\rm{t}}} = \frac{p}{{\gamma - 1}} + \frac{1}{2}\rho {u^2} + \frac{{{B^2}}}{{2{\mu _{\rm m}}}} , $

各向同性剪切应力张量为

$ {{\tau }} = - \mu \left[ {\nabla {{u}} + {{\left( {\nabla {{u}}} \right)}^{\rm{T}}}} \right] + \frac{2}{3}\mu \left( {\nabla \cdot {{u}}} \right){{I}} , $

各向异性剪切应力张量为

$ {{\tau }} = - {\mu _0}\nabla {{u}} - 3{\mu _{||}}\left( {{{bb}} - \frac{1}{3}{{I}}} \right)\left( {{{bb}} - \frac{1}{3}{{I}}} \right):\nabla {{u}} .$

上述方程中, $\rho$为流体密度, u为流体速度矢量, p为等离子体压力, I为二阶单位张量, B是磁感应强度, b为磁场单位矢量, $\mu_ {\rm{m}}$为磁导率, $\mu$为动力黏性系数(下标“0”表示各向同性分量, 下标“||”表示各向异性分量). 状态方程采用理想气体状态方程, 绝热指数$\gamma =1.4$.
本文使用有限体积法对方程进行离散, 采用PPM (piecewise parabolic method)三阶空间重构和HLLD (HLL-discontinuties)通量格式. 使用无方向分裂的CTU算法求解上述方程组, 应用CT算法处理磁场无散度约束.
为方便研究, 本文选择二维平面剪切层为研究对象, 如图1所示. 考虑两股在x方向上相向流动的自由剪切层流动, 假设流体是理想导电、具有各向异性黏性, 水平方向上流速分布为双曲正切间断,

$ {u_x}(y) = \frac{{\Delta U}}{2}\tanh \left(\frac{{2y}}{\delta }\right) , $

上下层速度差$ \Delta U = 0.5a $, 这里a为声速. 整个流场附加大小为B₀的流向均匀磁场,

${{{B}}_0} = {B_0}\hat x.$

而y方向(纵向)上附加高斯衰减的正弦速度扰动,

$\delta {u_y} = \varepsilon \exp \left( { - \frac{{{y^2}}}{{{\sigma ^2}}}} \right)\sin {k_x}x, $

其中$\varepsilon $为振幅, $\sigma $为波幅度衰减率, $\lambda$为波长, k_x =$ 2{\text{π}}/\lambda$是x方向上的扰动波数. 为防止扰动对流动干扰过大, 这里选择小振幅的扰动$\varepsilon = 0.03$, 并且在过渡层内快速衰减$\sigma = 2\delta $.
初始时刻, 来流中上、下层具有相同的初始压力p₀、密度$\rho_0$、比热比$\gamma = 1.4$、声速a = ($\gamma p_0/\rho_0$)^1/2, 以及阿尔文速度c_A = B₀/($\mu_ {\rm{m}}\rho_0$)^1/2. 为了获得涡结构最大分辨率, 设置扰动波长$\lambda = L$. 为了使得不稳定性增长较为明显, 选择相对较窄的过渡层$ \delta = L/50 $. 本文中, 采用L为单位长度, 声音跨越长度L的时间L/a为单位时间, ($\mu_ {\rm{m}}\rho_0$a²)^1/2作为单位磁场强度以无量纲化计算结果.
计算域尺寸为L_x × L_y = L × L, x方向两侧使用周期边界条件, y方向两侧采用出流边界条件. 网格采用均匀分布的笛卡尔网格, 经收敛测试后网格总数为1024 × 1024. 考虑到各向同性黏性对KH不稳定性有抑制作用, 因此选择雷诺数$Re=\rho_0\Delta UL/\mu$ = 10⁵、各向同性与各向异性雷诺数$Re_0=\rho_0\Delta UL/\mu_0$ = 10⁵, $Re_{||} =\rho_0\Delta UL/\mu_ {||} = $100分别进行数值模拟. 选择弱磁场情形(阿尔文马赫数$M_ {\rm{A}}= \Delta U/c_ {\rm{A}}= 10$)以便于观察黏性的影响. 为追踪上下层流体的演化, 在剪切层中加入示踪剂, 其浓度(质量分数)c由被动标量方程控制:

$\frac{{\partial \rho c}}{{\partial t}} + \nabla \cdot (\rho c{{u}}) = 0, $

并将其与MHD方程组耦合求解.
图 1 计算模型
Figure1. Sketch of the computational model.

图2为不同时刻剪切层示踪剂浓度场截图, 给出了各向同性黏性和各向异性黏性磁流体中KH不稳定性的形成与演变过程. 从图2可知, 各向同性黏性磁流体剪切层在扰动下界面发生卷起, 然而由于磁场的作用, 卷起结构被斜向拉长, 但不能继续卷起形成猫眼涡结构(图2(a)). 对于各向异性黏性情形(图2(b)), 剪切层仍然发生了卷起, 被斜向拉长, 但也不能形成典型的猫眼涡结构, 且卷起结构中心处出现了波纹状界面. 对比图2(a)和图2(b)可知, 图2(b)中卷起的圈数减少. 这些特征暗示各向异性黏性的加入具有致稳效果.
图 2 不同时刻自由剪切层示踪剂浓度场　(a)各向同性黏性(Re = 10⁵); (b) 各向异性黏性(Re₀ = 10⁵, Re_|| = 100)
Figure2. Tracer concentration field at different times: (a) Isotropic viscous case (Re = 10⁵); (b) anisotropic viscous case (Re₀ = 10⁵, Re_|| = 100).

由于黏性影响涡的扩散, 下面从涡结构演化的角度来分析其作用过程. 最近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)). 上述过程表明了黏性对剪切层涡结构的分裂增殖以及扩散作用, 使流动中出现了逆时针涡流, 从而感应出逆时针的速度场, 以对抗顺时针旋转的卷起.
图 3 不同时刻, 各向同性黏性(Re = 10⁵)算例涡结构云图及流线　(a) t = 5.5; (b) t = 7.0; (c) t = 7.5; (d) t = 8.0; (e) t = 9.0; (f) t = 10
Figure3. Rortex field with streamlines at different times in isotropic viscous fluid (Re = 10⁵): (a) t = 5.5; (b) t = 7.0; (c) t = 7.5; (d) t = 8.0; (e) t = 9.0; (f) t = 10.

图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)中心也形成了逆时针旋转的椭圆型长涡, 但上下两个涡流并未被强烈挤压.
图 4 不同时刻各向异性黏性(Re₀ = 10⁵, 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
Figure4. Rortex field with streamlines at different times in anisotropic viscous fluid (Re₀ =10⁵, 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为不同时刻各向同性与各向异性黏性算例中磁场强度和磁感线的分布. 先考察各向同性黏性情形(图5(a)). 在t = 5.5时, 剪切界面发生卷起, 由于磁冻结效应, 磁感线也随之卷起, 形成螺旋状结构. 卷起使磁感线被拉伸和弯曲, 导致磁场被放大, 因此图中能清晰地看到卷起的磁感线上磁场强度约为周围磁场的3倍以上. t = 8.0时, 界面进一步卷起并被斜向拉伸, 整体结构约45°倾斜, 图中可见磁感线卷起程度提高且在卷起外缘(红色区)磁场放大程度较强. 然而到t = 10时, 整个磁感线结构进一步被斜向拉长, 而且有向内被挤压的趋势, 因而磁场放大程度较强的(红色)区域分布较广. 但在卷起结构内部, 磁感线卷曲程度变得比t = 9.0时要低, 这反映了抗弯磁张力的稳定性作用.
图 5 不同时刻剪切层中磁场云图和磁感线　(a)各向同性黏性(Re = 10⁵); (b) 各向异性黏性(Re₀ = 10⁵, Re_|| = 100)
Figure5. Magnetic field with field lines at different times: (a) Isotropic viscous case (Re = 10⁵); (b) anisotropic viscous case (Re₀ = 10⁵, Re_|| = 100)

对于各向异性黏性情形(图5(b)), 可与各向同性黏性情形分别对比. 首先, 在结构方面, 各向异性黏性情形中磁感线卷起结构中的螺旋的圈数少于各向同性黏性情形, 中心处几乎没有卷曲形状. 其次, 在磁场强度方面, 各向异性黏性情形中卷起界面上磁场放大程度低于各向同性黏性情形, 且磁场被放大的区域分布范围不集中. 这表明, 各向异性黏性的加入, 使得沿磁感线方向上的剪切程度降低, 从而降低了磁感线卷曲程度, 导致磁场放大程度的下降, 使得流动变得稳定.
图6为流场中纵向总动能随时间的变化. 可见各向同性黏性算例和各向异性黏性算例都经历了线性增长与饱和后的非线性阶段. 然而不同的是, 各向异性黏性算例的线性增长率以及饱和水平都比各向同性黏性算例的低, 且线性增长时间稍微增加. 这从动能的角度表明, 各向异性黏性相比于各向同性黏性, 具有促进稳定的作用.
图 6 纵向总动能随时间的变化
Figure6. Evolution of the longitudinal total kinetic energy E_ky .

最近, Liu等^[27]研究发现, 磁场对KH不稳定性的致稳作用源于磁能局部放大所产生的横向磁压力(f_mp)和抗弯磁张力(f_mt), 下面分析这些磁场量的变化.
图7为平均磁能密度相对于初始时刻放大倍数随时间的变化. 可见, 在t = 3.5之前平均磁能密度几乎保持不变, 之后开始非线性上升, 直到大约在t = 9.5时第一次到达峰值. 其中各向同性黏性算例中平均磁能密度在峰值时放大到9倍左右, 而各向异性黏性算例中平均磁能密度则放大到6倍左右. 这正好解释了图4中局部磁场放大的区别.
图 7 平均磁能密度放大倍数随时间的变化
Figure7. Evolution of the amplification factor of average magnetic energy.

图8为平均横向磁压力和抗弯磁张力随时间的变化. 图8(a)中各向同性情形横向磁压力(红色实线)初期几乎线性上升到t = 4左右, 此后开始非线性增长, 直到t = 10左右饱和. 对于各向异性黏性算例(蓝色点线), 其曲线变化趋势与各向同性黏性算例相似, 但增长速度和幅值始终低于各向同性黏性情形. 这是因为各向异性黏性降低界面卷起程度, 减小磁能密度的放大程度, 从而降低了横向磁压力. 图8(b)中的抗弯磁张力则有所不同, 各向同性黏性与各向异性黏性算例对应曲线几乎重合, 二者都在t = 0—0.5左右快速线性上升, 然后非线性缓慢上升, 直到t = 10. 这表明各向异性黏性几乎对抗弯磁张力无影响.
图 8 平均横向磁压力和抗弯磁张力随时间的变化　(a) 横向磁压力; (b) 抗弯磁张力
Figure8. Evolution of the average transverse magnetic pressure and anti-bending magnetic tension: (a) Transverse magnetic pressure; (b) anti-bending magnetic tension.

图9为不同算例中流场总涡度拟能随时间的变化. 可见, 曲线在t = 0到t = 3.5左右几乎重合, 且都是缓慢下降. 随后, 各向同性黏性算例对应的曲线开始上升, 在t = 7之后几乎直线上升; 而各向异性黏性算例仍然沿之前趋势缓慢下降到t = 7.5左右才开始缓慢上升. 到t = 10时, 各向异性黏性算例中涡度拟能只有各向同性的25%左右. 这表明, 黏性和磁场在KH不稳定波的线性增长阶段发挥了稳定性作用, 缓慢降低了涡度拟能. 随后在非线性发展阶段, 由于黏性的扩散作用, 导致流动中生成了许多细小的涡, 而之前扰动动能则一部分转化为磁能, 另一部分则转化为涡度拟能. 各向异性黏性则由于其加强了磁感线方向的黏性, 使得沿着磁感线方向上的剪切速度变慢, 从而有效降低了涡量, 因此曲线上升速度不及各向同性黏性算例.
图 9 流场总涡度拟能随时间的变化
Figure9. Evolution of the total enstrophy.

利用CTU + CT算法对非理想磁流体方程组进行求解, 模拟了匀强平行磁场作用下, 具有各向异性黏性剪切层中的KH不稳定性, 研究了黏性各向异性对KH不稳定性的影响和作用机理.
模拟结果表明, 各向异性黏性具有稳定作用, 在同样情况下, 其致稳效果比各向同性黏性稍强. 其稳定作用是因为: 各向异性黏性中沿磁感线方向上的黏性系数比垂直于磁感线方向上的黏性系数要大得多, 导致界面卷起过程中沿着剪切面上的剪切速度变慢, 从而使卷起结构内部的卷曲圈数降低, 并且卷曲程度下降, 最终流动不再继续往不稳定方向发展.
在该过程中, 黏性使得卷起结构中小涡产生增殖, 特别是在边缘和内部分裂出数个逆时针的小涡. 这些小涡还会发生合并, 破坏了卷起结构常规增长, 从而使得流动变得稳定. 各向异性黏性对磁场的影响主要体现在降低横向磁压力方面, 对抗弯磁张力的影响较小.