磁场散度为零的变分重构法在日冕太阳风数值模拟中的应用

王浩鹏, 赵敬民, 吕佳锟, 柳晓静

1 中国科学院国家空间科学中心， 北京 100190 2 中国科学院大学地球与行星科学学院， 北京 100049

0 引言

太阳风是来自太阳的、持续的等离子体流(Abbo et al., 2016; Cranmer and Winebarger, 2019).太阳风充斥于日地空间并携带行星际磁场向远离太阳的方向运动(Li et al., 2021).剧烈的太阳爆发活动往往抛射出数以亿吨的高速日冕物质，朝向地球的高速运动的日冕物质能在几个小时到几天的时间内到达地球轨道，触发一系列灾害性空间天气.太阳风暴如日冕物质抛射引发的灾害性空间天气，对人类社会的诸多天基和地基高科技活动，如卫星导航、无线电通讯、宇航员太空活动、飞行器太空作业以及电力系统运转等造成巨大的经济损失.研究太阳活动机理，揭示太阳风暴传播演化规律，提前几小时到两三天时间预报出灾害空间天气，对维持人类的一系列高科技活动正常进行意义重大(冯学尚等, 2011,2013; Feng, 2020a).

日地空间尺度有上亿千米，太阳爆发活动在日地空间的传播演化包含了多时间和空间尺度的复杂物理过程，进行MHD数值建模成为目前唯一能够自洽的连接巨大的日地空间的研究手段.近几十年来，有一大批学者投身于用基于物理的三维MHD模型来进行日地空间的研究，极大地推动了太阳风MHD数值模拟在空间天气研究和应用领域中的发展(Feng 2020a).在这些模拟中，基于近似黎曼解的有限体积法应用比较广泛(Feng, 2020b).

在太阳附近的日冕区域磁场较强，等离子体β较小.在这些区域的等离子体团的磁压远大于热压，在数值模拟过程中，从能量导出的压强很容易出现非物理的负压.在这些小β区域，显式数值方法中的时间步长受柯朗稳定条件限制而取值较小(一般小于10-3h)，远远小于太阳风演化的时间尺度.因此，在设计用于模拟日冕太阳风结构的MHD模型时，应当充分考虑数值算法保持密度和压强为正数的性质(保正性)，计算效率的高效性.

HLL类型的黎曼算子具有较好的保正性.Harten-Lax-van Leer (HLL)黎曼算子最早由Harten等(1983)提出并用于求解双曲守恒系统，在HLL黎曼扇中解变量被两个快波分为两个外部状态和一个中间状态.近期，Wu和Shu(2019)在局部磁场散度为零的条件下严格证明了HLL黎曼算子的保正性，并用于二维磁流体力学问题中的数值模拟.由于HLL算子具有较强的稳定性而且计算量较少，本文采用具有保正性的HLL黎曼算子来计算通过网格界面的数值通量.

在计算效率方面，Feng等(2013)采用GPU高性能计算方式来提高计算速度，Wang等(2019)通过隐式时间积分方法扩大时间推进过程中的时间步长来提高收敛效率，进而显著减少了稳态背景太阳风数值模拟的计算量，将模拟时间从几天缩短到了几个小时甚至几十分钟.在Wang等(2019)和Feng等(2021)中分别用到了GMRES(Generalized Minimal Residual)算法和LU-SGS(lower-upper symmetric Gauss-Seidel)算法.GMRES最早由Saad和Schultz(1986)提出，是一种用于求解线性系统的全矩阵求解方法.相对于GMRES方法，另外一种用于求解线性系统的LU-SGS方法则是一种近似分解方法.LU-SGS方法将全矩阵[A]分解为上[L]、下[U]三角矩阵和对角阵[D]，并忽略了高阶无穷小项[L][D]-1[U]，即[A]=[L]+[D]+[U]≈([L]+[D])[D]-1([U]+[D])≡[L][D]-1[U].LU-SGS方法最早由Yoon和Jameson(1988)提出，具有无需求解全矩阵和占用存储空间少的优点，被广泛用于求解线性系统.本文将采用后向欧拉隐式时间积分法，并用LU-SGS方法求解每一次时间推进过程中的线性方程组.

此外，Feng等(2019)在最小二乘法的基础上加入平均磁场散度为零的限制设计出了保持全局磁场散度为零的磁场重构方法.相较于最小二乘法，变分重构方法(variable reconstruction method，VR)更加紧致，所需的重构模板更小(Liu et al., 2017; Li et al., 2018; Nishikawa, 2018; Wang et al., 2017)，有助于设计高阶高效率的重构方法.受此启发，本文基于VR方法，在重构磁场时加入平均磁场散度为零的限制条件，进而在变分重构法的基础上设计出了满足磁场散度为零的重构方法.

第1节介绍了本文日冕太阳风MHD模型中的控制方程.第2节介绍了本文中的网格系统和边界条件设置情况.第3节介绍了本文的数值方法，主要包括控制方程的离散以及本文提出的VR-GSP重构方法.第4节展示了本文基于VR-GSP方法的日冕太阳风MHD模型的模拟结果，并比较了分别基于本文VR-GSP方法和Feng等(2019)中的LSQ-GSP方法的MHD模型的模拟结果.第5节对本文中的方法做了总结，并对下一步的工作做了展望.

1 物理模型

日冕模型求解的是三维的理想MHD方程组.在近太阳表面的低等离子体β区域，由于磁压远大于热压，所以在由总能密度导出热压的浮点数计算中可能会得到非物理的负压值，从而导致模拟的失败.为避免负压值的出现，我们将磁场B分解为不随时间变化的势场B0=(B0x,B0y,B0z)T，以及随时间演化的小量B1=(B1x,B1y,B1z)T两部分，即B=B0+B1，进而得到如下形式的MHD方程组(Feng et al., 2010,2012a; Fuchs et al., 2010; Guo, 2015)：

(1)

其中，

(2)

(3)

2 网格系统和初边值条件

2.1 网格系统

图1 六片网格系统Fig.1 Six-component grid system

2.2 初边值条件

本文计算区域内的磁场B0由基于太阳全球振荡监测网(Global Oscillation Network Group GONG)观测到的光球层径向磁场数据(网址：https:∥gong2.nso.edu/archive/patch.pl?menutype=z)的势场模型生成，随时间变化的小量B1的初始量为0，初始密度ρ、径向速度vr和压强p由一维的Parker解(Parker, 1963)给出.太阳表面的质子数密度设置为1.5×108cm-3，太阳表面温度为1.3×106K.

内边界位于太阳表面，在此处施加回流边界条件(Groth et al., 2000; Feng et al., 2007, 2019)，根据当地的径向速度vr设置边界条件：

如果vr>0

如果vr≤0

外边界位于超声速、超阿尔芬速区域，本文直接将r2ρ，r2vr，vθ，rvφ，r2Br，Bθ，rBφ进行等值外推，其中vθ，vφ，Br，Bθ分别表示速度和磁场在θ-和φ-方向的分量.

3 数值方法

在本节中，我们主要介绍了控制方程的离散以及VR-GSP重构算法的推导.

3.1 空间积分

在网格单元celli上对方程(1)进行积分，得到方程的半离散格式(Feng, 2020b)：

(4)

其中，Γij表示网格单元celli和cellj的交界面，也表示交界面Γij的面积，nij表示界面Γij上由celli指向cellj的单位法向量，Fij·nij表示沿界面Γij法向量方向的对流通量.此时，方程(4)可简写为

(5)

由对流通量的旋转不变性(Tanaka, 1994)，得

(6)

在本文中，我们采用UnL和UnR分别表示T8ijUL和T8ijUR，下标L和下标R分别表示位于界面Γij左右两侧的状态.在上述表达式中，

是将(x,y,z)坐标系旋转到(n,t1,t2) 坐标系的旋转矩阵(Yalim, 2008; Feng, 2020b, and references therein)，其中，t1=(t1ij,x,t1ij,y,t1ij,z)T和t2=(t2ij,x,t2ij,y,t2ij,z)T表示界面Γij上两个正交的单位切向量，且与nij组成右手系的(n,t1,t2)直角坐标系.f(UnL,UnR)表示局部坐标系(n,t1,t2)中nij方向上的通量，因此，旋转不变性将多维问题转换成了局部的一维问题来进行计算.本文采用HLL黎曼算子来计算通量f(UnL,UnR)，即

SL和SR表示两个快波的速度，其具体取法请参考Wu和Shu(2019).

为计算方程(5)中残差项Ri的值，需要重构出网格celli的界面上的变量值.在本文中，我们采用二阶重构方法，从分别以网格celli和cellj为中心的模板中重构位于界面Γij的面心上的原始变量的值Wα=(ρα,uα,vα,wα,pα,B1x α,B1y α,B1z α)，α∈{L,R}.

3.2 线性最小二乘法重构

在网格celli中，ρ，u，v，w，p在位置x的重构表达式如下：

X∈{ρ，u，v，w，p}

3.3 VR-GSP重构方法

在以celli和cellj为中心的模板中，我们可以得到两网格共享界面Γij的面心xij处的磁场各分量值XL=Xi(xij)和XR=Xj(xij),X∈{Bx，By，Bz}，具体表达式如下：

(7)

(8)

其中，

Ω=diag(1,ωxLx,ωyLy,ωzLz)，

Lx，Ly，Lz表示网格长度尺度，ωx，ωy，ωz表示权重，{Vj}i表示仅包含网格单元celli和与celli共面的6个邻居网格的重构模板.利用中值定理计算方程(8)里的面积分，可得

(9)

其中，

ΔX=μjgj-μigi+(X|j-X|i)e1

由方程(9)得到

(10)

其中，

在方程(10)中，网格celli体心处的导数值gi和与celli共面的邻居网格cellj体心处的导数值gj均是待求解的未知量，因此计算区域所有网格中待求解的导数值是耦合在一起的.本文采用Gauss-Seidel迭代法(Wang et al., 2017)来求解这个耦合系统.

参考Feng等(2019)；Hopkins(2016)，网格celli中的磁场散度平均值可按如下方式计算

(11)

(12)

其中，

下标|i和|j分别表示相应变量在网格celli和cellj的体心处的取值.

(13)

其中，

然后通过迭代方式求解方程(13)，即

3.4 时间离散

在时间方向上，本文采用的是后向欧拉方法，将方程(5)写为

(14)

(15)

其中，Δhi表示网格celli的特征长度，本文取为celli的最小内切球半径，vn和cf,n分别表示celli的表面上沿外法向方向的等离子体速度与快磁声速，CFL表示柯朗数，在本文的隐式方法中取CFL=10.

从而得到如下线性系统

其中

D′(UnL,UnR)=D(UnL,UnR)+0.8λmax(UnR-UnL),

λmax=|vn|+cf,n.

进而得到下述线性系统，

(16)

4 模拟结果

在本节中，我们将本文提出的VR-GSP重构方法用于2234卡林顿周(CR 2234)的日冕太阳风模拟.为了说明该方法的有效性，我们将该方法和Feng等(2019)中的LSQ-GSP重构方法的模拟结果做了比较.为便于描述，我们将应用了VR-GSP和LSQ-GSP重构方法的隐式MHD模型分别简述为VR-GSP-MHD模型和LSQ-GSP-MHD模型.CR 2234 位于第24个太阳活动周的末期，从2020年8月11日开始，到2020年9月7日结束，本文选取GONG观测的2020年8月25日的光球层径向磁场作为模型的观测输入生成定态背景太阳风.图2展示的是CR 2234位于太阳表面的磁场强度等离子体β.在太阳表面大部分区域磁场强度大于1 Gs(高斯，1 Gs10-4T)，等离子体β值小于0.1.在局部区域，磁场强度最大值大于20 Gs，相应的等离子体β值小于3×10-4.在本文模拟中，时间步长Δt大约为0.0054 h，VR-GSP-MHD模型和LSQ-GSP-MHD模型分别经过6699次和7778次时间推进达到稳定状态，所需计算时间twall分别为6858 s和9152 s.模型计算均在国家超级计算天津中心的Th-1A超级计算机(https:∥www.nscc-tj.cn/)上进行，使用了12个计算节点，且每个计算节点配有两个Intel Xeon X5670 CPU(2.93 GHz, 6核).

图2 太阳表面的磁场强度(Gs)(a)和等离子体β(b)的概略图Fig.2 Synoptic maps of the magnetic field strength in units of Gauss (a) and the plasma β distribution (b) near the solar surface

4.1 太阳附近的日冕结构

根据汤姆逊散射原理得到的白光偏振图可以较好地反映近太阳区域的等离子体密度的分布(Hayes et al., 2001; Linker et al., 1999; Feng, 2020a)，图3分别展示了来自SOHO飞船(Brueckner et al., 1995)的广角光谱日冕仪观测结果(a,b)(下载网址 https:∥stereo-ssc.nascom.nasa.gov/browse/)和由VR-GSP-MHD模型模拟结果所得的白光偏振图(c, d),以及在所选子午面上的磁力线结构(e,f).在φ=180°-0°的子午面上，VR-GSP-MHD模型在φ=180°和φ=0°的半平面上分别重现了南纬10°和赤道附近沿太阳径向延伸的盔状明亮结构；在φ=270°和φ=90°的半平面上，分别重现了在赤道和北纬10°沿太阳径向方向延伸的明亮结构.(e)和(f)图中的磁力线结构进一步表明上述观测和模拟所得的明亮结构为冕流结构.

图3 SOHO飞船上的Large Angle and Spectrometric Coronagraph C2 (LASCO-C2)日冕仪观测的白光偏振亮度图(a,b)和使用VR-GSP-MHD模型模拟所得φ=180°-0°和φ=270°-90°两个所选子午面上的2.3～6Rs范围内的白光偏振亮度图(c,d) 以及在这两个子午面上1～5Rs范围内的磁力线分布图(e,f)Fig.3 White-light polarized brightness (pB) images observed from LASCO-C2/SOHO (a,b), and simulated images from 2.3 to 6Rs on the meridian planes of φ=180°-0° and φ=270°-90° that are synthesized from the results of the VR-GSP-MHD model (c,d), and the simulated magnetic field lines from 1 to 5Rs on the two selected meridian plans (e,f)

图4展示的是模拟得到的子午面上的太阳风速度和等离子体数密度的分布云图，云图上的带箭头实线表示磁力线，箭头表示磁力线的方向.图5展示了太阳风速度、等离子体数密度和温度在2.6Rs和20Rs的分布概略图.在图4中，低速高密度等离子体流主要集中在φ=180°半平面的以南纬10°为中心和φ=0°平半面的以北纬10°为中心，以及φ=270°半平面的以北纬10°为中心和φ=90°平面的以赤道为中心的大约30°范围的区域，高速低密度流则占据了中高纬区域，这体现出了极小年日冕太阳风的较为显著的结构特征.在图5中，磁场中性线在70°-210°经度范围有一个向南弯曲大约15°范围的圆弧结构，在其余经度范围则平缓分布，VR-GSP-MHD模型中的磁场中性线位置和低速高密度等离子体流比较集中的区域基本重合.该卡林顿周的磁中性线比较平缓，仅分布在赤道附近北纬10°和南纬10°的范围内，从而导致高速低密度流占据了中纬和高纬区域.在2.6Rs处，VR-GSP-MHD模型和PFSS模型产生的磁中性线位形基本一致，但存在细微的差异.由于PFSS模型将太阳表面和源表面间的磁场简化为势场，源表面外的磁场则根据源表面的径向磁场简单外推得到，并没有考虑磁场和等离子体流的相互作用，而VR-GSP-MHD模型所生成的定态太阳风中的磁场是初始的势场和等离子体流相互作用的结果，因此两者结果中的磁中性线位形存在一定的差异.

图4 使用VR-GSP-MHD模型模拟所得的在φ=180°-0°(a, c)和φ=270°-90°(b, d)两个子午面上的1～20Rs范围内径向速度vr(km·s-1) (a, b)，质子数密度(cm-3)以10为底的对数N(c, d).云图中黑色带有箭头的实线代表磁力线Fig.4 Magnetic field lines from 1 to 20Rs overlaid on contours of the radial speeds vr (km·s-1) (a, b) and the decadic logarithms of the proton number density (cm-3) (c, d) on the meridian planes of φ=180°-0° (a,c) and φ=270°-90° (b,d)

图5 由VR-GSP-MHD模型模拟得到的位于2.6Rs处(a,c,e)和20Rs处(b,d,f)的太阳风径向速度vr(km·s-1)(a,b)和质子数密度N(105cm-3)(c,d)以及等离子体温度T(105K)(e, f)的概略图.图中白色实线和黑色实线分别表示由VR-GSP-MHD模型和PFSS模型产生的磁中性线(Br=0)的位置.Fig.5 Synoptic maps of the radial speeds vr (km·s-1) (a,b), the proton number density (105cm-3) (c,d), and the temperature T(105K) (e,f) from the VR-GSP-MHD model at 2.6Rs (a,c,e) and 20Rs (b,d,f). The black solid and white solid lines denote the magnetic neutral lines (MNLs) from the PFSS and the VR-GSP-MHD models, respectively

图6展示了用VR-GSP-MHD和LSQ-GSP-MHD模型得到的径向速度和质子数密度的以10为底的对数在(θ，φ)=(5°，250°)的电流片附近和(θ，φ)=(-70°，250°)的冕洞区域沿太阳径向方向从1Rs到20Rs的分布图.两种方法模拟的分布曲线基本重合，从4Rs到20Rs，低速流从60 km·s-1增长到300 km·s-1，高速流从195 km·s-1增长到660 km·s-1，和由Pätzold等(1997)的日冕观测导出的结果一致.

图6 径向太阳风速度vr(km·s-1)(a)和质子数密度(cm-3)的以10为底的对数N(b)沿日心距离的分布.三角形和黑色实线分别表示由VR-GSP-MHD模型和LSQ-GSP-MHD模型得到的位于电流片附近的低速高密度流，菱形图标和黑色虚线分别表示由VR-GSP-MHD模型和LSQ-GSP-MHD模型得到的位于冕洞区域的高速低密度流Fig.6 Radial profiles of the radial speeds vr(km·s-1) (a) and the decadic logarithms of proton number density (cm-3) (b). The triangles and solid lines represent the profiles of the high-density and low-speed streams from VR-GSP-MHD and LSQ-GSP-MHD models, and the diamonds and dashed lines denote the variations of the low-density and high-speed streams from VR-GSP-MHD and LSQ-GSP-MHD models

图7 模拟和观测所得的径向速度vr(km·s-1)(a)和径向磁场极性(-1表示指向日心的方向，+1表示远离太阳的方向)(b)随日球经度的分布图.图中黑色实线表示Wind卫星在1 AU处的行星际观测值被映射到20Rs处的值，蓝色实线和红色虚线分别表示由VR-GSP-MHD模型和LSQ-GSP-MHD模型得到的20Rs处的径向速度和径向磁场极性Fig.7 Temporal profiles of the modeled and observed radial speeds vr(km·s-1) (a) and the radial magnetic field polarities (b) with “-1” denoting the field pointing to the Sun and "+1" the field directing away from the Sun at 20Rs. The black solid lines denote the mapped observations in situ interplanetary measurements, and the blue solid lines and red dashed lines denote the modeled results at 20Rs by VR-GSP-MHD and LSQ-GSP-MHD models, respectively

4.2 VR-GSP-MHD模型和LSQ-GSP-MHD模型性能的比较

参考Feng等(2019)，Hopkins(2016)和Powell等(1999)，我们分别定义单个网格单元celli和整个计算区域内的平均磁场散度误差Error(B)i和Error(B)ave为

图8和图9展示的均是关于磁场散度误差的结果.从图8中可以看出，在时间推进过程中，由VR-GSP和LSQ-GSP两种方法模拟得到的平均磁场散度误差保持在同一个量级上.起初，两种磁场散度消去方法中的平均磁场散度误差呈上升趋势，大约在10 h处达到峰值，分别为10-5和10-4.9，最终两种方法都能将磁场散度误差控制在10-5.5以下.图9显示的是稳态结果的磁场散度误差在子午面上的分布，显然，在大部分区域的磁场散度误差Error(B)i都小于10-7.

图8 由VR-GSP-MHD(虚线)和LSQ-GSP-MHD(实线)模型模拟的平均磁场散度误差以10为底的对数随松弛时间的演化图Fig.8 Temporal evolution profiles of Error(B)ave modeled by VR-GSP-MHD (dashed line) and LSQ-GSP-MHD (solid line) models

图9 由VR-GSP-MHD模型得到的相对磁场散度误差在所选取的两个子午面上1～20Rs范围内的分布云图Fig.9 Contour maps from 1 to 20Rs for the relative divergence errors of the magnetic field obtained from the VR-GSP-MHD model on the selected meridian planes

图10 由VR-GSP-MHD和LSQ-GSP-MHD模型得到质子数密度的相对差diffρ在φ=180°-0°(a)和φ=270°-90°(b)子午面1～20Rs范围内的分布云图Fig.10 Contour maps from 1 to 20Rs for the relative differences of proton number density diffρ on the meridional planes of φ=180°-0° (a) and φ=270°-90° (b)

图11 由VR-GSP-MHD和LSQ-GSP-MHD模型得到质子数密度的相对差在2.6Rs处(a)和20Rs处(b)的分布云图，白线表示VR-GSP-MHD模型模拟结果中的磁中性线Fig.11 Synoptic maps for the relative differences of the plasma density diffρ at 2.6Rs (a) and 20Rs (b). The white lines denote the MNLs from the VR-GSP-MHD model

在表1中，我们比较了在太阳风模拟中VR-GSP和LSQ-GSP的计算效率.表1分别列出了本文用VR-GSP-MHD和LSQ-GSP-MHD模型模拟CR 2234太阳风日冕结构达到稳定状态时所需的计算时间twall(单位：h)、迭代次数以及每次迭代所用的时间(单位：s).两者相比，我们可以看出VR-GSP不仅减少了每一步迭代计算所用的时间，也减少了模拟至稳定状态所需要的总的迭代次数，从而也说明了基于VR-GSP的MHD模拟的计算效率更高.

表1 VR-GSP-MHD和LSQ-GSP-MHD模型计算效率的对比Table 1 Comparison between computational efficiency of the VR-GSP-MHD and LSQ-GSP-MHD models

5 总结和展望

在本文中，我们基于变分重构提出一种可以消去磁场散度误差的重构方法，并将其简称为VR-GSP.通过上述对CR 2234期间的日冕太阳风的模拟表明，本文提出的VR-GSP方法不仅在消去磁场散度误差方面可以达到LSQ-GSP方法相近的效果，而且模拟出的日冕太阳风结构和观测结果基本一致.和LSQ-GSP方法相比，VR-GSP方法在一次时间推进中平均所需的计算时间减少了13.0%，达到稳态所需要的时间推进次数减少了13.9%，达到稳态所需计算时间减少了25.1%，可以在一定程度上提高计算效率.

由于变分重构法紧致，易于推广到高阶重构，因此本文中的方法有望进一步应用到高阶(>2阶)的三维太阳风MHD建模中.同时，把VR-GSP方法进一步用于数据驱动的太阳风数值模拟(Feng et al., 2012b,2015; 李会超等，2019)、行星际模拟(Guo et al., 2021)并且结合航天器联合观测(Zhang et al.,2022)进行MHD模型验证是我们正在考虑的下一步工作.

致谢本文工作使用了太阳全球振荡监测网(GONG)项目提供的数据.GONG项目由AURA, Inc. 运营的美国国家太阳观测台根据与美国国家科学基金会的合作协议管理.也使用了SOHO 飞船上的LASCO C2日冕仪和太阳动力学天文台(SDO)的大气成像组件 (AIA) 获得的数据.此外，我们还使用了NASA/GSFC 的OMNIWeb系统提供的数据，OMNI数据可由网址：https:∥cdaweb.gsfc.nasa.gov/index.html/下载.本文中的计算在国家超级计算天津中心的Th-1A超级计算机上完成.