耀斑大尺度电流片的湍流耗散和能谱分析∗

单 洁 叶 景 蔡强伟 林 隽

(1 中国科学院云南天文台昆明650011)(2 中国科学院大学北京100049)(3 洛阳师范学院空间物理研究所洛阳471934)(4 中国科学院天文大科学研究中心北京100012)

1 引言

磁重联是发生在磁化等离子体中的基本物理过程, 它广泛存在于实验室等离子体、地磁层、行星际空间和太阳活动等天体物理过程中. 磁重联是通过快速的磁场耗散, 将磁能转化为其他形式的能量(如等离子体的动能和热能以及高能带电粒子的动能), 并加热和加速等离子体, 同时引起等离子体中磁场拓扑结构发生变化的过程. 在太阳大气中,磁重联过程伴随一系列的太阳活动现象, 如日冕加热、太阳爆发和太阳风等[1−4]. 在地球磁尾中, 磁重联引起大尺度磁场位形改变, 引起剧烈的等离子体动力学过程, 造成磁暴和磁层亚暴[5]. 在超大质量黑洞与吸积盘系统中, 盘冕磁结构的爆发和随后的磁重联导致间歇性喷流与相应的吸积盘耀斑[6]; 而在恒星级黑洞和吸积盘系统中发生的前后相随的间歇性喷流会因为碰撞而发生磁重联, 这非常有可能是产生伽马暴的一种机制[7−8].在自转较慢, 而磁场超强(可达到1015Gs)的磁中子星(磁星)上, 有可能在磁星表面附近一些局部区域形成日珥那样的复杂磁场结构, 这些磁场结构失去平衡之后就会产生比太阳爆发要剧烈十几个量级的巨型罕见耀发(能量可超过1047erg). 在磁星表面的复杂结构的形成过程以及这些复杂结构失去平衡并产生巨型耀发的过程中, 磁重联都起着非同寻常的作用[9−11]. 磁重联时常也发生在小尺度的实验等离子体中[10−13].

经典的磁重联图像主要有两种: Sweet-Parker (SP)重联和Petschek重联. 对于发生在太阳爆发等剧烈活动中的能量释放过程来说, SP重联的速率太低, 不能解释驱动太阳爆发所需要的快速磁重联率. 这是因为SP模型中耗散区尺度太大, 日冕电阻率太小, 并且只有电阻耗散一种方式.

Petschek模型增加了新的耗散方式—慢模激波. 在这个模型中SP式的耗散区被限制在尺度远小于系统尺度的范围内, 耗散得到第1级提速; 同时由SP耗散引起的扰动以慢模激波的方式向外传播, 使得磁场的耗散能够更加迅速地进行, 耗散得到进一步的提速(Priest等[5]的文献对此有非常详细的讨论). 而近年来发展起来的湍流磁重联理论则进一步为增加耗散方式、加速磁重联率提供了一条新的路径[14−15]: SP式的电流片当中出现的等离子体不稳定性引起湍流, 形成众多的小尺度结构, 使得磁场的耗散速率在等离子体电阻率变化不大的情况下得到极大提高.

在具有有限电阻的电流片当中, 有3种宏观等离子体不稳定性可能发生: 重力模不稳定性、波纹模不稳定性以及撕裂模不稳定. 重力模是由于重力(或是等效重力)引起的密度分层所致, 波纹模则是由于电流片当中的温度不均匀引起的磁扩散不均匀所导致的,而撕裂模则能够在有电阻存在的长电流片中发生[12], 并且是3种不稳定性中最重要的一种. 由于它的长波特性, 撕裂模不稳定性引起的磁力线变形最平缓, 由此产生的恢复力最小, 因此它比其他两种不稳定性更容易发展. 其结果就是一整块的大尺度电流片被迅速“撕裂”成为一块包含有众多小尺度结构的“夹心饼干”. 在每个小尺度结构上同时发生的耗散要比在整个大尺度电流片上的耗散快得多, 这极大地增强了综合的耗散效果, 使得磁重联的整体速率要比简单的SP磁重联有了本质的飞跃. 这与使用并行计算的方法使得计算效率得到大幅提高类似, 最终的效果相当于在原来的SP电流片等离子体中增加了“超电阻(hyper-resistivity)”[14−16].

超电阻的出现源自于磁场被耗散的同时, 磁结构的尺度也被湍流有效地降低, 导致磁场耗散的非线性增强, 比大尺度区域中的耗散要更加有效[14−16]. 在这个过程中产生的随机分布(缠绕交织)的小尺度磁场造成了等效垂直于宏观平均场的动量传输以及反常的电子粘滞[17]. 在这样的情况下, Kolmogorov微观尺度(即湍流出现明显耗散的尺度)就不再单纯地由流体本身的特性所决定的, 而会与磁流体当中的局部发电机过程密切相关.

这些结果表明了大尺度磁重联过程的湍流性质, 揭示了磁流体在宏观尺度上发生耗散的物理本质, 也说明了湍流对加速磁重联过程的重要意义[18]. 另一方面, 在CME(Coronal Mass Ejection)-耀斑电流片中出现的湍流与传统的经典湍流图像又有些不同.首先, 在经典图像中, 湍流通过级联过程将动能或磁能从大尺度结构向小尺度结构转移,最终转化为无序的热能. 但是大尺度电流片中, 会出现两个或多个小尺度结构发生相互作用合并为大尺度结构的现象, 即发生了逆级联过程. 正反级联最后会达到一种动态平衡, 此时, 磁重联速率趋于常数, 而不再依赖于系统的初始状态及有关参数[19−21].

其次, 大尺度的电流片还可以容纳多种方式的耗散过程和结构同时出现, 并同时发挥作用. 比如, Mei等[22]的数值实验表明, 当电流片的长度足够长, 厚度足够大的时候,简单的SP结构、Petschek慢模激波以及多尺度的湍流结构都会出现在同一个电流片当中, 这在简单湍流导致的超电阻的基础上贡献了额外的耗散.

目前, 国际上已经有些研究组在等离子体物理的框架内发展完全3维的MHD模型方面取得了一些进展[23−24], 而且Daughton等[25]还在完全动力学的框架内对磁重联进行了3维模拟, 其结果表明电流片当中出现的许多3维结构在包含电阻的MHD模拟当中也会出现. Mei等[26]在完全3维的框架内研究了Titov-Demoulin磁场位形[27]的整体和局部演化, 也探讨了在这样的3维结构中发生磁重联的几何与物理特征. 他们发现系统失去平衡之后, 磁通量绳迅速向外运动并拉伸其周围的磁力线, 导致在通量绳之后形成一块电流片, 磁重联随后在其中发生[28]. 结果表明, 尽管这样形成的电流片是3维结构, 但它的整体结构基本上与2维的CME-耀斑电流片在第3个方向上的简单延伸没有本质的区别[22], 与2.5维磁重联电流片中的各种细节几乎完全一致[29]. 这说明片状结构是太阳爆发这类大尺度剧烈能量释放过程磁重联区域的普遍特征, 因此我们可以在2.5维甚至是2维的框架内详细研究这类磁重联过程.

在上述工作的基础上, 本工作将在不同磁雷诺数和空间分辨率的情况下研究湍流对磁重联率的影响, 考察湍流出现后, 磁重联电流片中耗散本质的变化, 探讨湍流能谱在不同条件下的表现. 本文的第2部分将对本工作所使用的模拟程序和方法进行简单介绍;第3部分给出计算结果, 并详细考察不同情况下磁重联率、湍流耗散、湍流能谱的表现和随时间的演化; 第4部分将讨论上述结果; 我们将在第5部分对本工作进行总结.

2 数值方法简介

本工作是Shen等[30]工作的延续. 在实验开始之前, 中线位于直角坐标系y轴上的电流片向上无限延伸, 分开了两边极性相反的、同样是向上无限延伸的、平行于电流片的磁场(见图1), 磁场与电流片的底端都位于太阳表面, 我们将这里设为坐标系的x轴, 电流片位于x= 0的位置上. 初始条件与Shen等[30]采用的初始条件相同, 系统在演化开始之前处于平衡状态. 制约系统演化的MHD方程与Shen等[30]使用的方程组完全一致, 边界条件也相同, 即底边界使用等离子体和磁场的系连条件, 而在其他3个方向上则使用开放条件.

在本工作中, 基本物理量的特征值取为: 磁场强度BN=50 Gs, 长度LN=1.4×108m, 数密度nN=1015m−3; 为了便于使用无量纲化, 因此进一步导出物理量特征值: 质量密度ρN= 1.67×10−12kg·m−3, Alfven速度vAN= 3.45×106m·s−1, 时标tN= 40.59 s,温度TN=7.21×108K. 由于计算资源的改善和计算技术的提高, 我们将考察表1所列的几种情况, 其中Ng是计算所用的格点数. 很明显本工作得到的计算结果的空间分辨率要明显优于Shen等[30]计算结果的空间分辨率.

图1 初始电流片和磁场位形, 其中颜色代表z方向上的电流密度JZ, 实线表示磁场线. 特征参数LN = 1.4×108 m、BN = 50 Gs、TN = 7.21×108 K.Fig.1 Initial current density and magnetic field, the color scale is current density along the z-axis JZ,solid lines are magnetic field lines. Characteristic parameter LN = 1.4×108 m, BN = 50 Gs,TN = 7.21×108 K.

表1 4组不同事例的相关参数Table 1 Four cases with different parameters

在本工作中, 我们使用ATHENA程序来进行模拟[31]. ATHENA是用于天体物理学中磁流体力学的基于格点计算的程序, 用96.1%的C语言、1.1%的Shell、1.0%的MATLAB (MATrix & LABoratory)、0.6%的M4、0.5%的Makefile、0.4%的IDL (Interactive Data Language)和0.3%的其他语言写成. 相比较之前的一些天体物理的计算程序, ATHENA有着更高的计算精度, 适合在计算时引入静态网格细化(SMR)和自适应网格细化(AMR), 用更低的计算资源获得更多更细致的结果. 目前ATHENA已经被广泛运用在天体物理和其他物理学科的数值模拟工作中. 在可压缩的绝热非粘性的理想磁流体力学环境中, ATHENA可以做1维、2维和3维的数值模拟工作. 通过修改或添加初始条件的参数和边界条件来模拟不同实际天体环境中的情况. 本工作中, ATHENA用来进行2维模拟, 而且使用静态网格来细化计算空间.

3 计算结果

在初始处于平衡状态的电流片上增加小扰动, 电流片两边的磁场就会因为互相吸引而靠近, 导致电流片被挤压而变薄. 在线性撕裂模不稳定性的框架内, 电流片的长宽比超过2π之后, 就会发生撕裂模不稳定性而使其内部结构产生变化[12]. 不过在实际情况中,线性撕裂模不稳定性不一定能够顺利发生, 电流片周围的非线性效应在一定程度上会抑制撕裂模不稳定性的发生. 比如, 在目前我们研究的系统中, 电流片底边界的系连条件就有可能阻碍线性撕裂模不稳定性的发展, 半无限长的尺度早就超过了线性撕裂模的临界长宽比. 即使在其他电流片有限长度的情况下, 长宽比达到50以上甚至100之后, 电流片中才有撕裂模不稳定性出现(参考文献[22]和[32]).

当电流片被挤压到临界状态, 磁岛开始在电流片最窄的部分出现, 原先电流片两边缓慢相向运动的等离子体和磁场突然加快了运动速度. 当越来越多的磁岛出现之后, 电流片中的磁重联进入大致稳定的快速进行阶段(图2).

3.1 磁重联率的变化

与其他的数值实验类似, 电流片最窄的部分变成了主X点, 在其两边出现的磁岛和磁重联出流分别向相反的两个方向运动; 同时在主X点附近有流体驻点出现, 在驻点周围, 流体的即时速度基本为零. 我们注意到驻点与主X点的空间位置交替变换. 由于主X点是磁重联外流和磁岛运动的分界点, 磁重联出流和磁岛都不会穿过主X点, 所以驻点在主X的哪一边, 磁岛就往哪边运动, 向不同方向运动的磁岛不会同时出现, 而是交替出现的. 这个景象再现了Shen等[30]的结果.

利用Shen等[30]和Mei等[22]的方法确定主X点的位置, 然后在其右侧沿着电流片的方向取宽为219 km、长为729 km的矩形区域, 通过计算区域中的入流速度和阿尔芬速度的平均值来估算磁重联率.

采用这个办法, 分别对不同分辨率和磁雷诺数的情况进行计算. 得出两种情况下磁重联率随时间的变化, 如图3所示. 可以看到, 在事例2中, 第1阶段的磁重联率很小, 一直保持在0.007以下(图3 (c)). 到时间t= 148左右, 磁重联率的值开始上升, 考察系统磁场结构和密度的演化过程, 我们将系统演化的密度图与其进行对比, 可以发现, 这个时间就是第1个磁岛开始出现的时间. 接下来磁重联率的值极快地增加, 在t= 148到t= 200这个时间段内, 磁重联率就从0.007攀升到了0.06. 接下来的演化时间里, 磁重联率随着磁岛的形成和演化开始在0.03和0.09之间波动, 与电流片中新磁岛出现的时间有关.

在事例3中(图3 (d)), 磁重联率变化的趋势和事例2大致相同, 只是第1阶段持续的时间更短. 第1阶段中, 即t= 62之前, 磁重联率小于0.005. 而随着第1个磁岛开始出现, 磁重联率攀升也略有加快, 在t= 62到t= 110这段时间, 从小于0.005攀升到0.075. 并且也在接下来的演化中随磁岛的产生和发展, 磁重联率在0.02到0.09之间上下波动.

图2 事例3中不同时刻的密度分布和磁场位形图, 其中颜色代表密度, 实线表示磁场线. 从(a)到(f)分别为t=28、80、150、309、440、524. 特征参数ρN = 1.67×10−12 kg·m−3、LN = 1.4×108 m、BN = 50 Gs、TN = 7.21×108 K.Fig.2 Initial density and magnetic field of different times in case 3 and the color represents density, solid lines are magnetic field lines. The time of (a) to (f) is t = 28,80,150,309,440,524, respectively.Characteristic parameter ρN = 1.67×10−12 kg·m−3, LN = 1.4×108 m, BN = 50 Gs, TN = 7.21×108 K.

事例1和事例4的表现与上述两种情况不太相同. 在事例1中, 通过考察系统演化图,我们发现第1个磁岛在t= 210时出现. 但重力影响改变了主X点及其周围磁场与等离子体结构的演化进程. 如图4 (a)所示, 事例1中的主X点位置t= 210到t= 443时段内发生了明显变化, 在这个时段内, 主X点上方有巨大的磁岛形成. 在重力的影响下磁岛带着主X点向下运动, 并最终落入下方的耀斑环中, 原先的主X点消失, 新的主X点瞬间转移到其他磁岛上方, 其位置也随着磁岛的运动而变化. 通过对比图3 (c)和图4 (c), 考虑到事例1的磁雷诺数Rm= 105, 事例2的磁雷诺数Rm= 106, 上述结果印证了Shen等[30]的结论, 即磁雷诺数越大, 磁重联过程越早进入湍流状态. 事例4与事例1类似. 图4中(c)、(d)所示, 在系统演化后期, 这两个情况的重联率都稳定在某个常数附近.

图3 事例2(左)和事例3(右)中的主X点位置和磁重联率随时间的分布图. (a)、(b)所示为主X点的位置; (c)、(d)所示为磁重联率MA随时间的变化.Fig.3 The position of principal X-point (PX-point) and MA in the case 2 (left) and case 3 (right). (a),(b) indicate the position of PX-point; (c), (d) show the MA with time.

最后我们将事例2和事例4进行对比, 这两种情况的磁雷诺数都是Rm= 106, 模拟网格数分别为3840×3840和7680×7680,事例4的细化程度是事例2的4倍. 因为事例4中磁岛受重力影响较早, 因此我们选取了0

3.2 湍流耗散分析

Mei等[22]将计算中获得的数据代入如下方程, 通过检验等式两边的平衡来评估计算当中数值耗散的大小:

其中A、v、B、η分别是磁矢势、流体速度、磁场强度和磁扩散率. 具体做法是在等式左边的耗散项中, 增加额外的耗散ηn, 于是(1)式变成:

移项、简化之后得到:

这里a=|∂tA −v×B+η∇×B|,b=|η∇×B|.

图4 事例1(左)和事例4(右)中的主X点位置和磁重联率随时间的分布图. (a)、(b)所示为主X点的位置; (c)、(d)所示为磁重联率MA随时间的变化.Fig.4 The position of PX-point and MA in the case 1 (left) and case 4 (right). (a), (b) indicate the position of PX-point; (c), (d) show the MA with time.

这里我们用“额外耗散”ηn来代表除经典Spitzer耗散之外的其他因素导致的耗散. 对于“额外耗散”, 我们首先想到的就是数值耗散. 在算法和程序给定之后, 数值耗散基本上就确定, 不会有大的变化, 当然不排除算法和程序本身存在缺陷, 导致原本不大的数值耗散被放大直至计算崩溃的情况. 根据我们这几年使用ATHENA程序的经验, 加上我们对一些算法的改善和补充, 在本工作的计算中还没有发生过这样的情况, 数值耗散保持在相对平稳的水平上.

除此之外, “额外耗散”还有可能包含了因为电流片中出现湍流之后产生的耗散. 很显然, (2)式描述的是非线性过程, 在湍流没有出现的时候, 非线性效应显现不出来, 当电流片受到等离子体不稳定性的影响, 出现湍流结构的时候, 这当中的非线性效应立刻就表现出来, 迅速放大电流片中的耗散效应, 相当于在原有的等离子体经典Spitzer电阻的基础上加了“超电阻”, 最终的效果就是磁重联速率显著增加(见文献[12–19,33,34]中的详细讨论). 因此, 我们可以利用(2)–(3)式开展两项工作, 首先当然是在电流片处于平流状态的时候, 估算数值耗散; 然后在电流片进入湍流状态之后, 考察其中的超电阻及其演化特征. 我们以主X点为中心, 在其周围选取边长为292 km的正方形区域, 提取区域内各点上的数据, 根据(3)式计算相应的ηn/η, 然后对这些结果取平均作为ηn/η的最终取值.

结果如图5和图6所示, 可以看到3组数值模拟中的额外耗散在系统演化初期都没有等于零的时候, 即数值误差在数值模拟试验中一定存在. 不同情况下额外耗散所占的比例也不同. 从图5中我们可以看见, 无论是事例2还是事例3, 系统演化的第1阶段额外耗散与经典耗散之比都很小, 事例2在0.05左右(图5 (a)), 事例3在5左右(图5 (b)). 在事例2中,t= 148左右, 即第1个磁岛开始出现时, 额外耗散与经典耗散之比也开始迅速增大, 增长幅度达30倍左右(图5 (a)). 接下来随着电流片中磁岛的生成和演化, 这两者的比值也出现上下波动. 这两个耗散之比在事例3当中的表现与事例2当中类似, 也是在湍流出现之前处于较小的水平, 随着湍流的出现而发生跃变, 增幅也在30倍左右(图5 (b)), 只是湍流起始的时刻、也就是额外耗散开始跃变的时刻明显早于事例2.

图5 额外耗散和物理耗散的比值大小随演化时间t的变化关系图. 其中(a)是事例2中的比值关系图, (a)中的插图为截取一段0 < t < 180; (b)是事例3中的比值关系图, (b)中的插图为截取一段0 < t < 120.Fig.5 The ratio of additional diffusion to physical diffusion with time. (a) is for case 2 with a illustration of 0 < t < 180 inside it; (b) is for case 3 with a illustration of 0 < t < 120 inside it.

这里的情况引起了我们的注意, 在事例3当中, 湍流开始前的额外耗散(此时以数值耗散为主)到达经典耗散的5倍左右. 这一方面说明数值耗散的不可避免, 另一方面也证实了我们早先的判断, 即数值耗散由具体的算法和程序决定, 一旦算法和程序确定之后,实际计算中的数值耗散一般会保持在比较平稳的水平上.

我们根据事例2和事例3的结果得出这个结论的依据如下: 在给定的磁化等离子体系统当中,Rm与系统的尺度和Alfven速度成正比, 与其中的磁扩散率(决定经典耗散率)成反比. 当尺度和Alfven速度不变时,Rm由磁扩散率η决定,η越大经典耗散越强,Rm越小;反之亦然. 事例2和事例3的基本磁场位形和强度都没有改变, 因此Rm的变化意味着经典耗散率的变化. 事例2和事例3当中的Rm相差100倍, 相当于其中的η相差100倍, 图5表明两个事例中的ηn/η正好相差100倍, 表明在两种情况下以数值耗散为主的ηn没有什么变化. 这也进一步说明ηn/η发生的跃变是由湍流引起, 相当于在磁重联过程中突然引入了“超级耗散”机制, 即超电阻.

为了比较相同磁雷诺数情况下不同空间分辨率对数值耗散的影响, 我们选取了事例2和事例4在湍流出现前的第1阶段50< t <150,来进行对比分析. 如图5 (a)中的插图和图6所示, 对ηn/η进行求平均值和求方差. 结果得出事例2中数值耗散和经典物理耗散比值的平均值为0.0970,而事例4中比值的平均值为0.0273, 小于事例2中的平均值, 事例2中的比值大小是事例4的3.55倍. 本次工作使用ATHENA 程序, 事例4的空间分辨率是事例2的2倍, 因此预测其数值耗散大小应该是事例2的1/4, 得到的结果与预期较为符合. 最后需要注意的是, 在本次工作中我们利用等时间间隔(∆t= 0.1tN, 其中∆t是输出数据的时间间隔)输出数据后验的方法来估计数值误差, 在一定程度上会高估真实的误差, 但提供了评估湍流耗散的新思路.

图6 事例4中, 在50 < t < 150的时间范围内, 额外耗散和物理耗散的比值大小随演化时间t的变化关系图.Fig.6 The ratio of additional diffusion to physical diffusion with time for a excerpt of 50 < t < 150 in case 4.

在结束这一节的工作之前我们需要指出, 对于使用的网格数(38402), 磁雷诺数对计算效果的影响在Rm= 106达到饱和, 同时从误差分析的图来看, 事例3的数值耗散是明显大于物理耗散的. 在本次工作中, 我们取Rm= 108的目的就是模拟由数值耗散主导的磁重联过程, 与由显式电阻主导的事例1和2进行对比. 在后面分析与磁雷诺数相关的结果时, 我们主要也是分析事例1、2和4之间的情形, 其间拿事例3的结果进行交叉对比. 结果表明在目前的数值计算环境中, 考察Rm=108的事例是没有意义的.

3.3 湍流能谱分析

自从撕裂模不稳定性[12]被发现40多年以来, 随着近20 yr来高性能计算能力的迅速增强, 由撕裂模不稳定性引起的磁流体动力学(MHD)湍流在加快磁场耗散的过程中所起的作用受到越来越多的关注[13,18,34−35]. B´arta等[20]曾利用高分辨率的数值模拟研究耀斑电流片中的分形结构和级联磁重联过程, 并给出沿着电流片方向的磁能分布的能谱指数为−2.14. Lazarian等[36]指出不是所有的混沌结构都是湍流, 湍流必然会导致能量级联到更小的尺度. 因此, 多段片状电流片和磁岛结构在2维模拟实验中可以被认为是湍流, 而磁岛合并过程中的反向级联不属于湍流. Dong等[37]也指出标准的惯性区域对应的谱指数为−1.5, 而当电流片中产生许多的磁岛结构以后谱指数会变得比Kolmogorov谱更陡峭. 因此, 研究磁能能谱有利于理解电流片动态演化时的结构和能量从大尺度到小尺度传递并最终耗散的特征.

我们沿着电流片中心线(x= 0)的方向, 取宽度293 km、长度2.8×104km的长条形区域, 对这个区域沿x方向求平均得到1维的磁能和动能分布, 通过傅里叶分析得到相应谱空间的能量分布特征. 针对磁重联发展的不同节点, 我们选取了两个代表性阶段: 撕裂模不稳定性开始发展阶段以及稳定重联的动态平衡阶段. 对这些时间点上的磁能和动能进行分析, 并得到相应的能谱随波数k的分布关系. 我们知道, 波数k是系统特征长度的倒数, 波数越小对应的结构尺度越大, 而波数越大则相应结构尺度越小. 因此, 能量随k的分布也就是不同尺度结构中的能量所占比例的分布. 本工作假设能谱分布是近似满足幂律分布的, 即E=a0k−γ, 其中E为能量密度,a0为常数. 将方程两边取对数, 得到lgE=lga0−γlgk, 对能谱的主体部分进行线性回归, 拟合出能谱分布图的斜率, 可以到谱指数γ. 图7–10是事例1–3的磁能和谱指数, 可以发现得到的斜率是负数, 相应的磁能谱指数γm和动能谱指数γk为正.

为了分析在同系统中不同重联阶段的能谱指数, 我们选取了事例3在t=80、309两个时刻的磁能分布曲线. 图7左列是在不同时间磁能沿x= 0的分布, 右列则是相应的傅里叶能谱. 在t= 80时, 撕裂模不稳定性开始发展, 磁重联进入非线性阶段, 其谱指数为γm= 2.31. 随着磁岛结构的逐步出现, 能量可以通过磁岛或者分形电流片结构高效地从大尺度向小尺度传递, 同时磁岛也可以通过融合达到能量从小尺度向大尺度的逆向传递, 与正向级联略有不同的是, 逆向级联过程中有磁重联发生, 因此逆向过程本身也会伴随着磁场的耗散. 在t= 309时, 重联处于正向级联和逆向级联的动态平衡, 同时, 随着磁岛的出现, 能谱的曲线出现了转折点, 在k= 421处谱指数发生偏转. 偏转前在较大尺度内, 即标准惯性区域, 对应的谱指数为γm= 1.51, 接近于经典湍流理论中的Kolmogorov谱(=5/3). 但由于众多磁岛的产生使得能谱在较小尺度上的分布比Kolmogorov谱更加陡峭, 其谱指数为γm= 2.17, 这个结果与Dong等[37]的结论基本一致. 值得注意的是, 线性回归得到的谱指数对于波数k的下边界选择十分敏感, 正如Clauset等[38]工作中讨论的. 我们进行的拟合区域大致相同, 主要集中在有明显幂律分布特征的大部分k值之间(一般大于10), 在这部分区域内能量分布较为密集, 并且由磁岛主导, 能够反映谱指数的大小.

在磁重联过程中, 等离子体的动能, 其中包括流体和磁岛的动能, 也会随着时间变化. Shen等[19]指出, 第1个磁岛的形成与能谱指数的变化不直接相关, 但与电流片中存在多个磁岛有关. 这是因为非线性磁重联过程的完全发展需要一定的时间. 级联磁重联在能谱指数的波动上会表现得很明显. 我们已经知道磁重联的演化过程分为两个阶段, 其分界的标志事件为第1个磁岛开始出现. 在第1个阶段中, 磁重联进行得十分缓慢,重联出流沿着电流片缓慢流动. 而当磁重联进入第2阶段后, 开始有撕裂模不稳定性出现, 形成磁岛并对磁重联过程产生极大的影响. 这个阶段中磁重联率会迅速增大. 因此, 对于不同时间点上动能的能谱进行分析也是十分有必要的. 我们也相应给出了事例3动能的分布曲线和能谱. 如图8所示, 动能能谱的分布相比于磁能能谱拟合的效果更好, 而且在k值较大的时候没有明显的偏转趋势. 在t=80、309时刻, 其谱指数γk分别为1.88和1.23.

图7 事例3中不同时刻磁能大小以及磁能能谱分布图. 其中图(a)、(b) t = 80,图(c)、(d) t = 309. 图(b)、(d)中实线画出了拟合图像, 谱指数的大小分别为2.31和2.17. 其中图(d)中k = 421为转折点, 转折点前的虚线所拟合出的谱指数为1.51.Fig.7 Magnetic energy and energy spectrum of different times in case 3. Panels (a) and (b) with the time of t = 80; Panels (c) and (d) with the time of t = 309. The solid lines in panels (b) and (d) are the fitting lines with indexes of 2.31 and 2.17, respectively. The turning point in panel (d) is k = 421, the dotted line is the fitting line with an index of 1.51 before the turning point.

接下来, 我们考察空间分辨率对能谱分布的影响, 图9分别给出了事例2在t= 613和事例4在t= 489的磁能和能谱分布. 如图9 (a)、(c)所示, 这两种情况下耀斑环高度基本相同, 且均进入了高速的稳定重联期. 在事例2中(图9 (b)), 磁能能谱在惯性区域内由谱指数γm= 1.68的过程占主导, 而在耗散区域由谱指数γm= 2.35的过程占主导. 两个区域的谱指数在k= 426处发生明显地偏转, 说明能量从大尺度向小尺度级联过程的惯性阶段在这里终止, 而耗散阶段开始. 通过计算发现, 此时发生耗散的尺度为657.3 km, 明显大于事例2单个网格的大小(72.9 km). 理论分析指出, 耗散一般只发生在动力学尺度,即粒子惯性半径附近(在太阳上, 通常只有几十米甚至几米的长度). 在我们的模型中, 典型电流片的耗散宽度约为154 km, 它可以通过测量横向穿过主X点的电流密度分布而得到. 又因为沿着电流片方向磁岛的长度往往是宽度的6倍左右, 所以磁能的分布理论上可以达到的耗散尺度约为6×154 =924 km. 模拟得到的耗散尺度与理论上相比在同等量级, 并且比理论小一些, 进一步验证了从这里耗散阶段开始的推论. 这些都说明磁能在电流片内的耗散可能发生在宏观的MHD尺度, 而不是必须得通过级联达到动力学尺度才能发生耗散. 另一方面, 通过只增加数值模拟网格数, 图(d)也给出了类似的结果.我们发现, 惯性区域的谱指数和耗散区域谱指数拐点发生在k= 481附近, 对应的尺度为582.1 km, 也远远大于单个网格的大小(36.5 km). 这表明, 不同的分辨率并没有明显改变Kolmogorov微尺度lko的大小.

在确定了系统中磁能进入耗散尺度不被空间分辨率所限制后, 我们接下来考察磁雷诺数对能谱分布的影响. 图10展示了事例1、2、3中t=449、337、268时的磁能能谱分布图, 这3种情况的空间分辨率一致, 磁雷诺数分别为Rm=105、Rm=106和Rm=108.由图10左列可见, 在这3个时间点上他们耀斑环的高度达到一致, 并且有相似的磁岛结构. 在图10右列, 我们发现3种情况惯性区域的磁能谱指数比较接近, 但拐点k的值略有区别. 在事例1中拐点k= 406, 意味着耗散开始发生的尺度约为688.0 km; 在事例2中拐点k= 427, 耗散开始发生的尺度约为655.7 km. 在经典的Kolmogorov湍流理论中, 磁流体当中的lko符合规律lko∼Rm−2/3. 因此, 事例2比事例1的耗散尺度更小,表现出一定的收敛趋势, 但由于数值误差的影响, 并没有严格地满足理论预期. 另外, 事例3相比事例2中磁雷诺数要高2个量级, 但拐点k= 464, 则耗散开始发生的尺度为603.4 km, 与事例2相差不大. 根据前面的误差分析得知, 事例3的数值耗散要比物理耗散(即Rm=108)大很多, 而事例2的数值耗散与物理耗散(即Rm=106)相当, 所以在实际计算中两者表现的真实磁雷诺数是接近的, 造成其耗散尺度也很接近.

4 讨论

我们的数值试验结果表明, 不同的磁雷诺数会影响撕裂模不稳定性发展的时间和湍动强度. 磁雷诺数越大磁岛出现的时间越早, 但湍流充分发展以后, 磁重联率就与初始设置的磁雷诺数无关. 实验发现, 在湍流出现前, 额外耗散与经典物理耗散之比保持在较小的水平, 随着磁岛出现而发生跃变. 这意味着湍流开始前的额外耗散以数值耗散为主, 且数值耗散会在接下来的实验中保持在比较平稳的水平上. 而出现磁岛后的湍流会使额外耗散与经典耗散的比值大幅度提高. 因此, 我们将在未来的工作中考察有多少能量是被湍流耗散掉的, 研究湍流在耀斑电流片内磁能转换所扮演的角色.

通过分析不同分辨率对磁能能谱拐点的影响, 我们发现耀斑电流片内的耗散尺度可能在宏观MHD尺度发生. 尤其是在电流片内存在慢模激波结构和磁岛相互融合的过程中会产生额外的耗散. 在以往的研究中, 惯性区域和耗散区域间的能谱上存在中间过渡区域, 能量可能会一边级联一边耗散. 这或许能为湍流磁重联过程在经典磁重联和经典湍流理论基础上提供全新的切入点. 而且, 耗散的尺度基本符合随磁雷诺数增大而减小的规律.

图8 事例3中不同时刻动能大小以及动能能谱分布图. 其中图(a)、(b) t = 80; 图(c)、(d) t = 309; 图(b)、(d)中斜率的拟合用实线标出, 谱指数分别为1.88、1.23.Fig.8 Kinetic energy and energy spectrum of different times in case 3. Panels (a) and (b) with the time of t = 80; Panels (c) and (d) with the time of t = 309; The solid lines in (b), (d) are the fitting lines with indexes of 1.88 and 1.23, respectively.

5 总结

在本次工作中, 我们利用2维非理想MHD数值模拟对耀斑电流片的形成和演化进行了分析. 共对4种情况进行了模拟计算, 分别是: 事例1, 磁雷诺数为Rm= 105, 格点数为3840×3840; 事例2, 磁雷诺数为Rm= 106, 格点数为3840×3840; 事例3, 磁雷诺数为Rm= 108, 格点数为3840×3840; 事例4, 磁雷诺数为Rm= 106, 格点数为7680×7680.通过不同的参数组合, 可以分别研究不同磁雷诺数和不同空间分辨率对系统演化的影响. 与Shen等[19,30]的工作相比, 增加了对影响湍流耗散和能谱分布的参数研究. 我们对模拟得到的结果进行分析和讨论, 主要讨论了3个部分的结果: 磁重联率的大小、湍流耗散以及湍流能谱的分析. 我们得到以下结论:

图9 事例2中t = 613和事例4中t = 489时的磁能大小以及磁能能谱分布图. 其中图(a)、(b)为事例2; 图(c)、(d)为事例4. 图(b)、(d)中转折点分别在k = 426和k = 481. 图(b)、(d)中的实线拟合出的转折后的谱指数分别为2.35和2.24.Fig.9 Magnetic energy and energy spectrum when t = 613 in case 2 and t = 489 in case 4. Panels (a)and (b) for case 2; Panels (c) and (d) for case 4. The turning points of panels (b) and (d) are k = 426 and k = 481. The solid lines in panels (b) and (d) are the fitting lines with indexes of 2.35 and 2.24 after the turning point, respectively.

(1)磁雷诺数和空间分辨率对模拟系统中磁重联率的大小影响不大, 但会对初始磁岛出现的时间产生影响. 磁雷诺数越大, 系统就可以越快产生初始磁岛, 越早进入非线性阶段. 当系统的磁重联从线性阶段进入非线性阶段, 湍流的发展会使磁重联率产生明显的抬升, 但最终稳定的重联率大小与初始条件关联性不大;

(2)我们用“额外耗散”来代表经典Spitzer耗散之外的、由其他因素导致的耗散, 主要包括数值耗散和出现湍流后产生的耗散. 数值耗散由具体的算法和程序决定, 一旦这两者确定后, 数值耗散一般会保持在比较平稳的水平上. 磁雷诺数相同时, 不同的空间分辨率会对数值耗散产生一定的影响, 分辨率越高这个影响就越小. 额外耗散与经典耗散之比ηn/η发生的跃变是由湍流引起的. 当系统演化至非线性阶段, 出现湍流结构的时候, 电流片中的非线性效应由此表现出来, 迅速放大电流片中的耗散效应, 并最终表现为磁重联率显著增加. 不同磁雷诺数对额外耗散的水平影响不大;

图10 事例1、2、3中t=449、337、268时的磁能大小和能谱分布图, 其中图(a)、(b)为事例1, 图(c)、(d)为事例2,图(e)、(f)为事例3. 图(b)、(d)、(f)中的转折点分别为k=406、427、464, 转折后的能谱指数分别为2.29、2.89和2.85.Fig.10 Magnetic energy and energy spectrum when t = 449 in case 1, t = 337 in case 2, t = 268 in case 3. Panels (a) and (b) for case 1; Panels (c) and (d) for case 2; Panels (e) and (f) for case 3. The turning points in panels (b), (d) and (f) are corresponding to k = 406,427,464, and the solid lines in panels (b),(d), and (f) are the fitting lines with indexes of 2.29, 2.89 and 2.85 after the turning point, respectively.

(3)在重联处于线性阶段时, 大尺度磁场结构占主导, 谱指数相对较小; 当撕裂模不稳定性开始发展, 重联进入非线性阶段, 谱指数随着磁结构的成长上升至较高的水平; 而随着磁岛结构大量地产生, 能量稳定高效地从大尺度向小尺度传递, 同时磁岛融合等过程使能量从小尺度向大尺度逆向传递, 此时的谱指数绝对值降低, 磁能分布较为平缓.对动能能谱的拟合相比于磁能能谱拟合的效果更好一些, 但谱指数相对较小, 意味着动能在各个尺度的分布较均匀;

(4)当系统进入到高速稳定重联期后, 能谱分布曲线表明磁能能谱的谱指数会发生明显的偏转, 能量从惯性区域进入到耗散区域. 通过对4个事例中能谱的拐点进行分析,我们发现磁能在电流片中的耗散可能会发生在宏观MHD尺度, 不同的分辨率对磁能耗散的尺度影响不大, 改变网格划分精度不会改变系统进入耗散尺度的大小. 基于同一种网格划分, 磁雷诺数较高时对应的耗散尺度更小, 呈现一定的收敛趋势, 基本符合经典的湍流理论.

致谢此项工作得到中国科学院云南天文台计算太阳物理实验室支持, 相关数值计算在实验室计算平台上完成.