磁場(chǎng)散度為零的變分重構(gòu)法在日冕太陽(yáng)風(fēng)數(shù)值模擬中的應(yīng)用

王浩鵬, 趙敬民, 呂佳錕, 柳曉靜

1 中國(guó)科學(xué)院國(guó)家空間科學(xué)中心， 北京 100190 2 中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049

0 引言

太陽(yáng)風(fēng)是來(lái)自太陽(yáng)的、持續(xù)的等離子體流(Abbo et al., 2016; Cranmer and Winebarger, 2019).太陽(yáng)風(fēng)充斥于日地空間并攜帶行星際磁場(chǎng)向遠(yuǎn)離太陽(yáng)的方向運(yùn)動(dòng)(Li et al., 2021).劇烈的太陽(yáng)爆發(fā)活動(dòng)往往拋射出數(shù)以億噸的高速日冕物質(zhì)，朝向地球的高速運(yùn)動(dòng)的日冕物質(zhì)能在幾個(gè)小時(shí)到幾天的時(shí)間內(nèi)到達(dá)地球軌道，觸發(fā)一系列災(zāi)害性空間天氣.太陽(yáng)風(fēng)暴如日冕物質(zhì)拋射引發(fā)的災(zāi)害性空間天氣，對(duì)人類社會(huì)的諸多天基和地基高科技活動(dòng)，如衛(wèi)星導(dǎo)航、無(wú)線電通訊、宇航員太空活動(dòng)、飛行器太空作業(yè)以及電力系統(tǒng)運(yùn)轉(zhuǎn)等造成巨大的經(jīng)濟(jì)損失.研究太陽(yáng)活動(dòng)機(jī)理，揭示太陽(yáng)風(fēng)暴傳播演化規(guī)律，提前幾小時(shí)到兩三天時(shí)間預(yù)報(bào)出災(zāi)害空間天氣，對(duì)維持人類的一系列高科技活動(dòng)正常進(jìn)行意義重大(馮學(xué)尚等, 2011,2013; Feng, 2020a).

日地空間尺度有上億千米，太陽(yáng)爆發(fā)活動(dòng)在日地空間的傳播演化包含了多時(shí)間和空間尺度的復(fù)雜物理過程，進(jìn)行MHD數(shù)值建模成為目前唯一能夠自洽的連接巨大的日地空間的研究手段.近幾十年來(lái)，有一大批學(xué)者投身于用基于物理的三維MHD模型來(lái)進(jìn)行日地空間的研究，極大地推動(dòng)了太陽(yáng)風(fēng)MHD數(shù)值模擬在空間天氣研究和應(yīng)用領(lǐng)域中的發(fā)展(Feng 2020a).在這些模擬中，基于近似黎曼解的有限體積法應(yīng)用比較廣泛(Feng, 2020b).

在太陽(yáng)附近的日冕區(qū)域磁場(chǎng)較強(qiáng)，等離子體β較小.在這些區(qū)域的等離子體團(tuán)的磁壓遠(yuǎn)大于熱壓，在數(shù)值模擬過程中，從能量導(dǎo)出的壓強(qiáng)很容易出現(xiàn)非物理的負(fù)壓.在這些小β區(qū)域，顯式數(shù)值方法中的時(shí)間步長(zhǎng)受柯朗穩(wěn)定條件限制而取值較小(一般小于10-3h)，遠(yuǎn)遠(yuǎn)小于太陽(yáng)風(fēng)演化的時(shí)間尺度.因此，在設(shè)計(jì)用于模擬日冕太陽(yáng)風(fēng)結(jié)構(gòu)的MHD模型時(shí)，應(yīng)當(dāng)充分考慮數(shù)值算法保持密度和壓強(qiáng)為正數(shù)的性質(zhì)(保正性)，計(jì)算效率的高效性.

HLL類型的黎曼算子具有較好的保正性.Harten-Lax-van Leer (HLL)黎曼算子最早由Harten等(1983)提出并用于求解雙曲守恒系統(tǒng)，在HLL黎曼扇中解變量被兩個(gè)快波分為兩個(gè)外部狀態(tài)和一個(gè)中間狀態(tài).近期，Wu和Shu(2019)在局部磁場(chǎng)散度為零的條件下嚴(yán)格證明了HLL黎曼算子的保正性，并用于二維磁流體力學(xué)問題中的數(shù)值模擬.由于HLL算子具有較強(qiáng)的穩(wěn)定性而且計(jì)算量較少，本文采用具有保正性的HLL黎曼算子來(lái)計(jì)算通過網(wǎng)格界面的數(shù)值通量.

在計(jì)算效率方面，F(xiàn)eng等(2013)采用GPU高性能計(jì)算方式來(lái)提高計(jì)算速度，Wang等(2019)通過隱式時(shí)間積分方法擴(kuò)大時(shí)間推進(jìn)過程中的時(shí)間步長(zhǎng)來(lái)提高收斂效率，進(jìn)而顯著減少了穩(wěn)態(tài)背景太陽(yáng)風(fēng)數(shù)值模擬的計(jì)算量，將模擬時(shí)間從幾天縮短到了幾個(gè)小時(shí)甚至幾十分鐘.在Wang等(2019)和Feng等(2021)中分別用到了GMRES(Generalized Minimal Residual)算法和LU-SGS(lower-upper symmetric Gauss-Seidel)算法.GMRES最早由Saad和Schultz(1986)提出，是一種用于求解線性系統(tǒng)的全矩陣求解方法.相對(duì)于GMRES方法，另外一種用于求解線性系統(tǒng)的LU-SGS方法則是一種近似分解方法.LU-SGS方法將全矩陣[A]分解為上[L]、下[U]三角矩陣和對(duì)角陣[D]，并忽略了高階無(wú)窮小項(xiàng)[L][D]-1[U]，即[A]=[L]+[D]+[U]≈([L]+[D])[D]-1([U]+[D])≡[L][D]-1[U].LU-SGS方法最早由Yoon和Jameson(1988)提出，具有無(wú)需求解全矩陣和占用存儲(chǔ)空間少的優(yōu)點(diǎn)，被廣泛用于求解線性系統(tǒng).本文將采用后向歐拉隱式時(shí)間積分法，并用LU-SGS方法求解每一次時(shí)間推進(jìn)過程中的線性方程組.

此外，F(xiàn)eng等(2019)在最小二乘法的基礎(chǔ)上加入平均磁場(chǎng)散度為零的限制設(shè)計(jì)出了保持全局磁場(chǎng)散度為零的磁場(chǎng)重構(gòu)方法.相較于最小二乘法，變分重構(gòu)方法(variable reconstruction method，VR)更加緊致，所需的重構(gòu)模板更小(Liu et al., 2017; Li et al., 2018; Nishikawa, 2018; Wang et al., 2017)，有助于設(shè)計(jì)高階高效率的重構(gòu)方法.受此啟發(fā)，本文基于VR方法，在重構(gòu)磁場(chǎng)時(shí)加入平均磁場(chǎng)散度為零的限制條件，進(jìn)而在變分重構(gòu)法的基礎(chǔ)上設(shè)計(jì)出了滿足磁場(chǎng)散度為零的重構(gòu)方法.

第1節(jié)介紹了本文日冕太陽(yáng)風(fēng)MHD模型中的控制方程.第2節(jié)介紹了本文中的網(wǎng)格系統(tǒng)和邊界條件設(shè)置情況.第3節(jié)介紹了本文的數(shù)值方法，主要包括控制方程的離散以及本文提出的VR-GSP重構(gòu)方法.第4節(jié)展示了本文基于VR-GSP方法的日冕太陽(yáng)風(fēng)MHD模型的模擬結(jié)果，并比較了分別基于本文VR-GSP方法和Feng等(2019)中的LSQ-GSP方法的MHD模型的模擬結(jié)果.第5節(jié)對(duì)本文中的方法做了總結(jié)，并對(duì)下一步的工作做了展望.

1 物理模型

日冕模型求解的是三維的理想MHD方程組.在近太陽(yáng)表面的低等離子體β區(qū)域，由于磁壓遠(yuǎn)大于熱壓，所以在由總能密度導(dǎo)出熱壓的浮點(diǎn)數(shù)計(jì)算中可能會(huì)得到非物理的負(fù)壓值，從而導(dǎo)致模擬的失敗.為避免負(fù)壓值的出現(xiàn)，我們將磁場(chǎng)B分解為不隨時(shí)間變化的勢(shì)場(chǎng)B0=(B0x,B0y,B0z)T，以及隨時(shí)間演化的小量B1=(B1x,B1y,B1z)T兩部分，即B=B0+B1，進(jìn)而得到如下形式的MHD方程組(Feng et al., 2010,2012a; Fuchs et al., 2010; Guo, 2015)：

(1)

其中，

(2)

(3)

2 網(wǎng)格系統(tǒng)和初邊值條件

2.1 網(wǎng)格系統(tǒng)

圖1 六片網(wǎng)格系統(tǒng)Fig.1 Six-component grid system

2.2 初邊值條件

本文計(jì)算區(qū)域內(nèi)的磁場(chǎng)B0由基于太陽(yáng)全球振蕩監(jiān)測(cè)網(wǎng)(Global Oscillation Network Group GONG)觀測(cè)到的光球?qū)訌较虼艌?chǎng)數(shù)據(jù)(網(wǎng)址：https:∥gong2.nso.edu/archive/patch.pl?menutype=z)的勢(shì)場(chǎng)模型生成，隨時(shí)間變化的小量B1的初始量為0，初始密度ρ、徑向速度vr和壓強(qiáng)p由一維的Parker解(Parker, 1963)給出.太陽(yáng)表面的質(zhì)子數(shù)密度設(shè)置為1.5×108cm-3，太陽(yáng)表面溫度為1.3×106K.

內(nèi)邊界位于太陽(yáng)表面，在此處施加回流邊界條件(Groth et al., 2000; Feng et al., 2007, 2019)，根據(jù)當(dāng)?shù)氐膹较蛩俣葀r設(shè)置邊界條件：

如果vr>0

如果vr≤0

外邊界位于超聲速、超阿爾芬速區(qū)域，本文直接將r2ρ，r2vr，vθ，rvφ，r2Br，Bθ，rBφ進(jìn)行等值外推，其中vθ，vφ，Br，Bθ分別表示速度和磁場(chǎng)在θ-和φ-方向的分量.

3 數(shù)值方法

在本節(jié)中，我們主要介紹了控制方程的離散以及VR-GSP重構(gòu)算法的推導(dǎo).

3.1 空間積分

在網(wǎng)格單元celli上對(duì)方程(1)進(jìn)行積分，得到方程的半離散格式(Feng, 2020b)：

(4)

其中，Γij表示網(wǎng)格單元celli和cellj的交界面，也表示交界面Γij的面積，nij表示界面Γij上由celli指向cellj的單位法向量，F(xiàn)ij·nij表示沿界面Γij法向量方向的對(duì)流通量.此時(shí)，方程(4)可簡(jiǎn)寫為

(5)

由對(duì)流通量的旋轉(zhuǎn)不變性(Tanaka, 1994)，得

(6)

在本文中，我們采用UnL和UnR分別表示T8ijUL和T8ijUR，下標(biāo)L和下標(biāo)R分別表示位于界面Γij左右兩側(cè)的狀態(tài).在上述表達(dá)式中，

是將(x,y,z)坐標(biāo)系旋轉(zhuǎn)到(n,t1,t2) 坐標(biāo)系的旋轉(zhuǎn)矩陣(Yalim, 2008; Feng, 2020b, and references therein)，其中，t1=(t1ij,x,t1ij,y,t1ij,z)T和t2=(t2ij,x,t2ij,y,t2ij,z)T表示界面Γij上兩個(gè)正交的單位切向量，且與nij組成右手系的(n,t1,t2)直角坐標(biāo)系.f(UnL,UnR)表示局部坐標(biāo)系(n,t1,t2)中nij方向上的通量，因此，旋轉(zhuǎn)不變性將多維問題轉(zhuǎn)換成了局部的一維問題來(lái)進(jìn)行計(jì)算.本文采用HLL黎曼算子來(lái)計(jì)算通量f(UnL,UnR)，即

SL和SR表示兩個(gè)快波的速度，其具體取法請(qǐng)參考Wu和Shu(2019).

為計(jì)算方程(5)中殘差項(xiàng)Ri的值，需要重構(gòu)出網(wǎng)格celli的界面上的變量值.在本文中，我們采用二階重構(gòu)方法，從分別以網(wǎng)格celli和cellj為中心的模板中重構(gòu)位于界面Γij的面心上的原始變量的值Wα=(ρα,uα,vα,wα,pα,B1x α,B1y α,B1z α)，α∈{L,R}.

3.2 線性最小二乘法重構(gòu)

在網(wǎng)格celli中，ρ，u，v，w，p在位置x的重構(gòu)表達(dá)式如下：

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

3.3 VR-GSP重構(gòu)方法

在以celli和cellj為中心的模板中，我們可以得到兩網(wǎng)格共享界面Γij的面心xij處的磁場(chǎng)各分量值XL=Xi(xij)和XR=Xj(xij),X∈{Bx，By，Bz}，具體表達(dá)式如下：

(7)

(8)

其中，

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

Lx，Ly，Lz表示網(wǎng)格長(zhǎng)度尺度，ωx，ωy，ωz表示權(quán)重，{Vj}i表示僅包含網(wǎng)格單元celli和與celli共面的6個(gè)鄰居網(wǎng)格的重構(gòu)模板.利用中值定理計(jì)算方程(8)里的面積分，可得

(9)

其中，

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

由方程(9)得到

(10)

其中，

在方程(10)中，網(wǎng)格celli體心處的導(dǎo)數(shù)值gi和與celli共面的鄰居網(wǎng)格cellj體心處的導(dǎo)數(shù)值gj均是待求解的未知量，因此計(jì)算區(qū)域所有網(wǎng)格中待求解的導(dǎo)數(shù)值是耦合在一起的.本文采用Gauss-Seidel迭代法(Wang et al., 2017)來(lái)求解這個(gè)耦合系統(tǒng).

參考Feng等(2019)；Hopkins(2016)，網(wǎng)格celli中的磁場(chǎng)散度平均值可按如下方式計(jì)算

(11)

(12)

其中，

下標(biāo)|i和|j分別表示相應(yīng)變量在網(wǎng)格celli和cellj的體心處的取值.

(13)

其中，

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

3.4 時(shí)間離散

在時(shí)間方向上，本文采用的是后向歐拉方法，將方程(5)寫為

(14)

(15)

其中，Δhi表示網(wǎng)格celli的特征長(zhǎng)度，本文取為celli的最小內(nèi)切球半徑，vn和cf,n分別表示celli的表面上沿外法向方向的等離子體速度與快磁聲速，CFL表示柯朗數(shù)，在本文的隱式方法中取CFL=10.

從而得到如下線性系統(tǒng)

其中

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

λmax=|vn|+cf,n.

進(jìn)而得到下述線性系統(tǒng)，

(16)

4 模擬結(jié)果

在本節(jié)中，我們將本文提出的VR-GSP重構(gòu)方法用于2234卡林頓周(CR 2234)的日冕太陽(yáng)風(fēng)模擬.為了說(shuō)明該方法的有效性，我們將該方法和Feng等(2019)中的LSQ-GSP重構(gòu)方法的模擬結(jié)果做了比較.為便于描述，我們將應(yīng)用了VR-GSP和LSQ-GSP重構(gòu)方法的隱式MHD模型分別簡(jiǎn)述為VR-GSP-MHD模型和LSQ-GSP-MHD模型.CR 2234 位于第24個(gè)太陽(yáng)活動(dòng)周的末期，從2020年8月11日開始，到2020年9月7日結(jié)束，本文選取GONG觀測(cè)的2020年8月25日的光球?qū)訌较虼艌?chǎng)作為模型的觀測(cè)輸入生成定態(tài)背景太陽(yáng)風(fēng).圖2展示的是CR 2234位于太陽(yáng)表面的磁場(chǎng)強(qiáng)度等離子體β.在太陽(yáng)表面大部分區(qū)域磁場(chǎng)強(qiáng)度大于1 Gs(高斯，1 Gs10-4T)，等離子體β值小于0.1.在局部區(qū)域，磁場(chǎng)強(qiáng)度最大值大于20 Gs，相應(yīng)的等離子體β值小于3×10-4.在本文模擬中，時(shí)間步長(zhǎng)Δt大約為0.0054 h，VR-GSP-MHD模型和LSQ-GSP-MHD模型分別經(jīng)過6699次和7778次時(shí)間推進(jìn)達(dá)到穩(wěn)定狀態(tài)，所需計(jì)算時(shí)間twall分別為6858 s和9152 s.模型計(jì)算均在國(guó)家超級(jí)計(jì)算天津中心的Th-1A超級(jí)計(jì)算機(jī)(https:∥www.nscc-tj.cn/)上進(jìn)行，使用了12個(gè)計(jì)算節(jié)點(diǎn)，且每個(gè)計(jì)算節(jié)點(diǎn)配有兩個(gè)Intel Xeon X5670 CPU(2.93 GHz, 6核).

圖2 太陽(yáng)表面的磁場(chǎng)強(qiáng)度(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 太陽(yáng)附近的日冕結(jié)構(gòu)

根據(jù)湯姆遜散射原理得到的白光偏振圖可以較好地反映近太陽(yáng)區(qū)域的等離子體密度的分布(Hayes et al., 2001; Linker et al., 1999; Feng, 2020a)，圖3分別展示了來(lái)自SOHO飛船(Brueckner et al., 1995)的廣角光譜日冕儀觀測(cè)結(jié)果(a,b)(下載網(wǎng)址 https:∥stereo-ssc.nascom.nasa.gov/browse/)和由VR-GSP-MHD模型模擬結(jié)果所得的白光偏振圖(c, d),以及在所選子午面上的磁力線結(jié)構(gòu)(e,f).在φ=180°-0°的子午面上，VR-GSP-MHD模型在φ=180°和φ=0°的半平面上分別重現(xiàn)了南緯10°和赤道附近沿太陽(yáng)徑向延伸的盔狀明亮結(jié)構(gòu)；在φ=270°和φ=90°的半平面上，分別重現(xiàn)了在赤道和北緯10°沿太陽(yáng)徑向方向延伸的明亮結(jié)構(gòu).(e)和(f)圖中的磁力線結(jié)構(gòu)進(jìn)一步表明上述觀測(cè)和模擬所得的明亮結(jié)構(gòu)為冕流結(jié)構(gòu).

圖3 SOHO飛船上的Large Angle and Spectrometric Coronagraph C2 (LASCO-C2)日冕儀觀測(cè)的白光偏振亮度圖(a,b)和使用VR-GSP-MHD模型模擬所得φ=180°-0°和φ=270°-90°兩個(gè)所選子午面上的2.3～6Rs范圍內(nèi)的白光偏振亮度圖(c,d) 以及在這兩個(gè)子午面上1～5Rs范圍內(nèi)的磁力線分布圖(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展示的是模擬得到的子午面上的太陽(yáng)風(fēng)速度和等離子體數(shù)密度的分布云圖，云圖上的帶箭頭實(shí)線表示磁力線，箭頭表示磁力線的方向.圖5展示了太陽(yáng)風(fēng)速度、等離子體數(shù)密度和溫度在2.6Rs和20Rs的分布概略圖.在圖4中，低速高密度等離子體流主要集中在φ=180°半平面的以南緯10°為中心和φ=0°平半面的以北緯10°為中心，以及φ=270°半平面的以北緯10°為中心和φ=90°平面的以赤道為中心的大約30°范圍的區(qū)域，高速低密度流則占據(jù)了中高緯區(qū)域，這體現(xiàn)出了極小年日冕太陽(yáng)風(fēng)的較為顯著的結(jié)構(gòu)特征.在圖5中，磁場(chǎng)中性線在70°-210°經(jīng)度范圍有一個(gè)向南彎曲大約15°范圍的圓弧結(jié)構(gòu)，在其余經(jīng)度范圍則平緩分布，VR-GSP-MHD模型中的磁場(chǎng)中性線位置和低速高密度等離子體流比較集中的區(qū)域基本重合.該卡林頓周的磁中性線比較平緩，僅分布在赤道附近北緯10°和南緯10°的范圍內(nèi)，從而導(dǎo)致高速低密度流占據(jù)了中緯和高緯區(qū)域.在2.6Rs處，VR-GSP-MHD模型和PFSS模型產(chǎn)生的磁中性線位形基本一致，但存在細(xì)微的差異.由于PFSS模型將太陽(yáng)表面和源表面間的磁場(chǎng)簡(jiǎn)化為勢(shì)場(chǎng)，源表面外的磁場(chǎng)則根據(jù)源表面的徑向磁場(chǎng)簡(jiǎn)單外推得到，并沒有考慮磁場(chǎng)和等離子體流的相互作用，而VR-GSP-MHD模型所生成的定態(tài)太陽(yáng)風(fēng)中的磁場(chǎng)是初始的勢(shì)場(chǎng)和等離子體流相互作用的結(jié)果，因此兩者結(jié)果中的磁中性線位形存在一定的差異.

圖4 使用VR-GSP-MHD模型模擬所得的在φ=180°-0°(a, c)和φ=270°-90°(b, d)兩個(gè)子午面上的1～20Rs范圍內(nèi)徑向速度vr(km·s-1) (a, b)，質(zhì)子數(shù)密度(cm-3)以10為底的對(duì)數(shù)N(c, d).云圖中黑色帶有箭頭的實(shí)線代表磁力線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)的太陽(yáng)風(fēng)徑向速度vr(km·s-1)(a,b)和質(zhì)子數(shù)密度N(105cm-3)(c,d)以及等離子體溫度T(105K)(e, f)的概略圖.圖中白色實(shí)線和黑色實(shí)線分別表示由VR-GSP-MHD模型和PFSS模型產(chǎn)生的磁中性線(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模型得到的徑向速度和質(zhì)子數(shù)密度的以10為底的對(duì)數(shù)在(θ，φ)=(5°，250°)的電流片附近和(θ，φ)=(-70°，250°)的冕洞區(qū)域沿太陽(yáng)徑向方向從1Rs到20Rs的分布圖.兩種方法模擬的分布曲線基本重合，從4Rs到20Rs，低速流從60 km·s-1增長(zhǎng)到300 km·s-1，高速流從195 km·s-1增長(zhǎng)到660 km·s-1，和由P?tzold等(1997)的日冕觀測(cè)導(dǎo)出的結(jié)果一致.

圖6 徑向太陽(yáng)風(fēng)速度vr(km·s-1)(a)和質(zhì)子數(shù)密度(cm-3)的以10為底的對(duì)數(shù)N(b)沿日心距離的分布.三角形和黑色實(shí)線分別表示由VR-GSP-MHD模型和LSQ-GSP-MHD模型得到的位于電流片附近的低速高密度流，菱形圖標(biāo)和黑色虛線分別表示由VR-GSP-MHD模型和LSQ-GSP-MHD模型得到的位于冕洞區(qū)域的高速低密度流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 模擬和觀測(cè)所得的徑向速度vr(km·s-1)(a)和徑向磁場(chǎng)極性(-1表示指向日心的方向，+1表示遠(yuǎn)離太陽(yáng)的方向)(b)隨日球經(jīng)度的分布圖.圖中黑色實(shí)線表示W(wǎng)ind衛(wèi)星在1 AU處的行星際觀測(cè)值被映射到20Rs處的值，藍(lán)色實(shí)線和紅色虛線分別表示由VR-GSP-MHD模型和LSQ-GSP-MHD模型得到的20Rs處的徑向速度和徑向磁場(chǎng)極性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)，我們分別定義單個(gè)網(wǎng)格單元celli和整個(gè)計(jì)算區(qū)域內(nèi)的平均磁場(chǎng)散度誤差Error(B)i和Error(B)ave為

圖8和圖9展示的均是關(guān)于磁場(chǎng)散度誤差的結(jié)果.從圖8中可以看出，在時(shí)間推進(jìn)過程中，由VR-GSP和LSQ-GSP兩種方法模擬得到的平均磁場(chǎng)散度誤差保持在同一個(gè)量級(jí)上.起初，兩種磁場(chǎng)散度消去方法中的平均磁場(chǎng)散度誤差呈上升趨勢(shì)，大約在10 h處達(dá)到峰值，分別為10-5和10-4.9，最終兩種方法都能將磁場(chǎng)散度誤差控制在10-5.5以下.圖9顯示的是穩(wěn)態(tài)結(jié)果的磁場(chǎng)散度誤差在子午面上的分布，顯然，在大部分區(qū)域的磁場(chǎng)散度誤差Error(B)i都小于10-7.

圖8 由VR-GSP-MHD(虛線)和LSQ-GSP-MHD(實(shí)線)模型模擬的平均磁場(chǎng)散度誤差以10為底的對(duì)數(shù)隨松弛時(shí)間的演化圖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模型得到的相對(duì)磁場(chǎng)散度誤差在所選取的兩個(gè)子午面上1～20Rs范圍內(nèi)的分布云圖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模型得到質(zhì)子數(shù)密度的相對(duì)差diffρ在φ=180°-0°(a)和φ=270°-90°(b)子午面1～20Rs范圍內(nèi)的分布云圖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模型得到質(zhì)子數(shù)密度的相對(duì)差在2.6Rs處(a)和20Rs處(b)的分布云圖，白線表示VR-GSP-MHD模型模擬結(jié)果中的磁中性線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中，我們比較了在太陽(yáng)風(fēng)模擬中VR-GSP和LSQ-GSP的計(jì)算效率.表1分別列出了本文用VR-GSP-MHD和LSQ-GSP-MHD模型模擬CR 2234太陽(yáng)風(fēng)日冕結(jié)構(gòu)達(dá)到穩(wěn)定狀態(tài)時(shí)所需的計(jì)算時(shí)間twall(單位：h)、迭代次數(shù)以及每次迭代所用的時(shí)間(單位：s).兩者相比，我們可以看出VR-GSP不僅減少了每一步迭代計(jì)算所用的時(shí)間，也減少了模擬至穩(wěn)定狀態(tài)所需要的總的迭代次數(shù)，從而也說(shuō)明了基于VR-GSP的MHD模擬的計(jì)算效率更高.

表1 VR-GSP-MHD和LSQ-GSP-MHD模型計(jì)算效率的對(duì)比Table 1 Comparison between computational efficiency of the VR-GSP-MHD and LSQ-GSP-MHD models

5 總結(jié)和展望

在本文中，我們基于變分重構(gòu)提出一種可以消去磁場(chǎng)散度誤差的重構(gòu)方法，并將其簡(jiǎn)稱為VR-GSP.通過上述對(duì)CR 2234期間的日冕太陽(yáng)風(fēng)的模擬表明，本文提出的VR-GSP方法不僅在消去磁場(chǎng)散度誤差方面可以達(dá)到LSQ-GSP方法相近的效果，而且模擬出的日冕太陽(yáng)風(fēng)結(jié)構(gòu)和觀測(cè)結(jié)果基本一致.和LSQ-GSP方法相比，VR-GSP方法在一次時(shí)間推進(jìn)中平均所需的計(jì)算時(shí)間減少了13.0%，達(dá)到穩(wěn)態(tài)所需要的時(shí)間推進(jìn)次數(shù)減少了13.9%，達(dá)到穩(wěn)態(tài)所需計(jì)算時(shí)間減少了25.1%，可以在一定程度上提高計(jì)算效率.

由于變分重構(gòu)法緊致，易于推廣到高階重構(gòu)，因此本文中的方法有望進(jìn)一步應(yīng)用到高階(>2階)的三維太陽(yáng)風(fēng)MHD建模中.同時(shí)，把VR-GSP方法進(jìn)一步用于數(shù)據(jù)驅(qū)動(dòng)的太陽(yáng)風(fēng)數(shù)值模擬(Feng et al., 2012b,2015; 李會(huì)超等，2019)、行星際模擬(Guo et al., 2021)并且結(jié)合航天器聯(lián)合觀測(cè)(Zhang et al.,2022)進(jìn)行MHD模型驗(yàn)證是我們正在考慮的下一步工作.

致謝本文工作使用了太陽(yáng)全球振蕩監(jiān)測(cè)網(wǎng)(GONG)項(xiàng)目提供的數(shù)據(jù).GONG項(xiàng)目由AURA, Inc. 運(yùn)營(yíng)的美國(guó)國(guó)家太陽(yáng)觀測(cè)臺(tái)根據(jù)與美國(guó)國(guó)家科學(xué)基金會(huì)的合作協(xié)議管理.也使用了SOHO 飛船上的LASCO C2日冕儀和太陽(yáng)動(dòng)力學(xué)天文臺(tái)(SDO)的大氣成像組件 (AIA) 獲得的數(shù)據(jù).此外，我們還使用了NASA/GSFC 的OMNIWeb系統(tǒng)提供的數(shù)據(jù)，OMNI數(shù)據(jù)可由網(wǎng)址：https:∥cdaweb.gsfc.nasa.gov/index.html/下載.本文中的計(jì)算在國(guó)家超級(jí)計(jì)算天津中心的Th-1A超級(jí)計(jì)算機(jī)上完成.