外日球?qū)蛹げㄊ录囊痪S磁流體力學(xué)數(shù)值模擬

郭孝城，周昱成，王 赤，李 暉

（1.中國(guó)科學(xué)院 國(guó)家空間科學(xué)中心 空間天氣學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2.中國(guó)科學(xué)院大學(xué) 地球和行星科學(xué)學(xué)院，北京 100049）

引 言

源于太陽(yáng)表面的超聲速太陽(yáng)風(fēng)等離子體向外徑向膨脹時(shí)會(huì)與包含星際等離子體和中性原子的鄰近星際介質(zhì)相互作用形成一泡狀結(jié)構(gòu)即日球?qū)?，它們的交界面就是日球?qū)禹擺1]。在經(jīng)歷漫長(zhǎng)的旅程后，“旅行者1號(hào)”（ Voyager 1）和“旅行者2號(hào)”（ Voyager 2）飛船分別在距太陽(yáng)94 AU和84 AU左右穿越終止激波，最終在120 AU左右穿越日球?qū)禹斶M(jìn)入星際空間[2-3]。太陽(yáng)風(fēng)事件，包括共轉(zhuǎn)太陽(yáng)風(fēng)作用區(qū)（Corotating Interaction Regions，CIRs）[4]和行星際日冕拋射事件（Interplanetary Coronal Mass Ejections，ICMEs）[5]，在日球?qū)觾?nèi)傳輸時(shí)會(huì)引起系列反應(yīng)，比如飛船發(fā)現(xiàn)在幾十AU處太陽(yáng)風(fēng)事件傳輸時(shí)會(huì)逐漸匯聚合并形成更大尺度的等離子體脈沖或激波事件[6-8]

“旅行者1號(hào)”在2012年8月穿越日球?qū)禹敽?，磁?qiáng)計(jì)在2 012.92（十進(jìn)制年份）和2014.66年（十進(jìn)制年份）分別探測(cè)到兩個(gè)前向激波，在2013.35年探測(cè)到一個(gè)疑似后向激波[9]。這些激波被認(rèn)為與“旅行者1號(hào)”的電場(chǎng)波動(dòng)儀器探測(cè)到的射電爆信號(hào)相關(guān)[10]。

學(xué)者嘗試對(duì)這些星際空間激波的來(lái)源做分析。Liu等[11]首先采用了一維磁流體力學(xué)數(shù)值模擬研究了太陽(yáng)風(fēng)事件在外日球?qū)訁^(qū)的傳播和演化，提出“旅行者1號(hào)”在2013.35年觀測(cè)到的射電爆和激波信號(hào)可能源自Wind飛船在2012年3月探測(cè)到的一系列的行星際日冕物質(zhì)拋射事件。然而，他們的模擬沒(méi)有考慮星際等離子體的入流，所以結(jié)果沒(méi)有終止激波和日球?qū)忧蕝^(qū)的形成，太陽(yáng)風(fēng)始終為超聲速；在計(jì)算激波傳播時(shí)類(lèi)比地球弓激波上下游的情況并對(duì)預(yù)定終止激波以上區(qū)域的激波速度進(jìn)行人為校正。Fermo等[12]采用更復(fù)雜的全球磁流體力學(xué)數(shù)值模擬研究了太陽(yáng)風(fēng)事件在外日球?qū)拥膫鞑ズ脱莼诰嗳? AU處引入OMNI（https://omniweb.gsfc.nasa.gov）的太陽(yáng)風(fēng)觀測(cè)數(shù)據(jù)，發(fā)現(xiàn)產(chǎn)生的日球?qū)禹敽图げㄎ恢帽葘?shí)際觀測(cè)到的數(shù)據(jù)要遠(yuǎn)約30 AU。后來(lái)，Kim等[13]擴(kuò)展了該工作，在1 AU的低緯度地區(qū)采用OMNI數(shù)據(jù)，但高緯度區(qū)采用經(jīng)驗(yàn)給出的高速太陽(yáng)風(fēng)數(shù)值[14]。結(jié)果顯示多重的太陽(yáng)風(fēng)共轉(zhuǎn)作用區(qū)對(duì)第一個(gè)前向激波的形成具有重要的作用，而第二個(gè)前向激波可能與太陽(yáng)風(fēng)共轉(zhuǎn)作用區(qū)和行星際日冕拋射事件的共同作用有關(guān)。

當(dāng)“旅行者1號(hào)”在星際空間探測(cè)到這3個(gè)激波時(shí)，“旅行者2號(hào)”此時(shí)在內(nèi)日球?qū)忧蕝^(qū)的太陽(yáng)風(fēng)里。Richandson等[15]通過(guò)“旅行者2號(hào)”的觀測(cè)數(shù)據(jù)發(fā)現(xiàn)鞘區(qū)內(nèi)的太陽(yáng)風(fēng)大尺度壓力脈沖與星際空間的激波具有一定的關(guān)聯(lián)性，經(jīng)過(guò)對(duì)波速的計(jì)算認(rèn)為5個(gè)合并太陽(yáng)風(fēng)作用區(qū)可能最終驅(qū)動(dòng)了相應(yīng)激波在星際空間的形成。然而，由于“旅行者1號(hào)”和“旅行者2號(hào)”并不沿著相同的方向飛行，這些壓力脈沖事件與星際激波的因果關(guān)聯(lián)仍有待確認(rèn)。本文將通過(guò)一個(gè)一維磁流體力學(xué)數(shù)值模擬對(duì)這個(gè)因果關(guān)系進(jìn)行分析和討論。

1 數(shù)值模型

首先假設(shè)太陽(yáng)風(fēng)在赤道面附近為沿著徑向方向運(yùn)動(dòng)的球?qū)ΨQ(chēng)流，這里忽略了太陽(yáng)風(fēng)在傳播過(guò)程中可能受到的來(lái)自側(cè)面方向太陽(yáng)風(fēng)的影響。此時(shí)理想磁流體力學(xué)方程組在球坐標(biāo)系下可以簡(jiǎn)化為沿著徑向 r方向的一維問(wèn)題為

其中：ρ、u、B和E分別是太陽(yáng)風(fēng)等離子體的密度、速度、磁場(chǎng)強(qiáng)度和能量密度。

總壓 pT和能量密度 E的具體形式為

其中：p是等離子熱壓。

考慮了太陽(yáng)風(fēng)等離子體離子與星際中性原子的電荷交換，方程組右邊的 QN、 QMr、 QM?和 QE分別表示因電荷交換引起的質(zhì)量、徑向動(dòng)量、方位角動(dòng)量和能量的變化[16-17]。太陽(yáng)的重力效應(yīng)也被考慮，其中G 是萬(wàn)有引力常數(shù)，Ms是太陽(yáng)質(zhì)量。采用MUSCL數(shù)值格式和有限體積方法數(shù)值求解上述方程組[18]，并運(yùn)用拓展型HLLC黎曼算子完成對(duì)方程通量的計(jì)算[19]，在時(shí)間推進(jìn)上采用龍格-庫(kù)塔格式，最終數(shù)值格式在空間和時(shí)間上均具有二級(jí)精度。

數(shù)值模擬區(qū)域范圍為r方向從1～200 AU，共1萬(wàn)個(gè)網(wǎng)格點(diǎn)。采用非均勻網(wǎng)格，設(shè)N為總格點(diǎn)數(shù)，則格面值ri(i=3,···,N ?1)由代數(shù)式為

其余四個(gè)格面值 r1、r2、rN、rN+1由插值給出；而格心值為兩相鄰格面值的中點(diǎn)；δ 參數(shù)表征最大網(wǎng)格與最小網(wǎng)格的差異程度，δ →0時(shí)網(wǎng)格趨向均勻網(wǎng)格，這里取 δ = 2.0。這種網(wǎng)格規(guī)模能夠使得模擬獲得較高的分辨率，比如在內(nèi)邊界網(wǎng)格尺寸約0.006 2 AU，在終止激波（～90 AU）附近約0.024 AU，在日球?qū)禹敚ā?20 AU）約為0.03 AU。初始時(shí)刻，內(nèi)邊界處太陽(yáng)風(fēng)取值為：數(shù)密度5 / cm3，速度ur=400 km/s，溫度1.5×105K，磁場(chǎng)強(qiáng)度 Br=2.8 nT。星際介質(zhì)方面，這里只考慮中性原子，其分布形式為[20]

其中：nH0=0.15/cm3為臨近星際空間的中性原子密度初始值； r0=4 AU為星際中性原子進(jìn)入日球?qū)拥纳疃?，用距離日心的距離表征。

簡(jiǎn)化計(jì)，模擬沒(méi)有考慮星際等離子體，因此太陽(yáng)風(fēng)與星際介質(zhì)相互作用的形式主要以電荷交換形式出現(xiàn)。

從初態(tài)開(kāi)始，太陽(yáng)風(fēng)經(jīng)過(guò)5年的演化達(dá)到準(zhǔn)穩(wěn)態(tài)，如圖1的點(diǎn)線(xiàn)所示。注意到，由于沒(méi)有星際等離子體的作用，終止激波和日球?qū)禹斁鶝](méi)有出現(xiàn)。由于太陽(yáng)風(fēng)離子與中性原子的電荷交換作用導(dǎo)致太陽(yáng)風(fēng)被加熱，在外日球?qū)訁^(qū)域其溫度沿著徑向并沒(méi)有出現(xiàn)絕熱冷卻的現(xiàn)象，而是在5 AU以上逐漸升高[21]。進(jìn)一步以此為初態(tài)，在1 AU的內(nèi)邊界處代之以2010—2017年間的OMNI、STEREO A和B（以下簡(jiǎn)寫(xiě)STA和STB）等飛船的觀測(cè)數(shù)據(jù)。因?yàn)镾TA在2014年9月—2015年11月間部分?jǐn)?shù)據(jù)缺失，而STB自2014年10月起已經(jīng)失去聯(lián)系，所以缺失的數(shù)據(jù)這里以上述恒定的初態(tài)條件代替。后面討論的外日球?qū)蛹げㄊ录?015年之前，假設(shè)激波速度500 km/s，那么從1～120 AU以上區(qū)域激波傳播需要1.14 a左右，因此可以認(rèn)為這種代替的負(fù)面影響是小的。圖1中的實(shí)線(xiàn)表示以O(shè)MNI數(shù)據(jù)為輸入，模擬數(shù)據(jù)在2013.42年時(shí)刻的物理量徑向分布，可見(jiàn)在太陽(yáng)風(fēng)事件沿著徑向方向向外傳播，出現(xiàn)的激波信號(hào)將會(huì)在120 AU以上的星際空間被探測(cè)到。

圖1 沿著徑向方向的太陽(yáng)風(fēng)參數(shù)Fig.1 The profiles of solar wind plasma along radial direction

2 模擬結(jié)果

需要說(shuō)明的是，在模擬中太陽(yáng)風(fēng)被簡(jiǎn)化為接近赤道面的球?qū)ΨQ(chēng)流?！奥眯姓?號(hào)”在2010—2017年間在HGI（Heliopsheric Inertial Coordinate）坐標(biāo)系下的緯度變化為34.4°N～ 34.8°N，“旅行者2號(hào)”從32.3°N～28.8°N；經(jīng)度變化上，“旅行者1號(hào)”是173.6°S～175.1°S，“旅行者2號(hào)”是216.8°S～218.3°S；可見(jiàn)兩艘飛船均離赤道面有一定的距離。但由于太陽(yáng)風(fēng)結(jié)構(gòu)在外日球?qū)訁^(qū)域傳播存在大區(qū)域傳播的特性，在經(jīng)度和緯度上覆蓋面廣[22]。這里假設(shè)模擬的結(jié)果能夠一定程度上反映兩艘飛船觀測(cè)得到的等離子體大尺度結(jié)構(gòu)的特征。

2.1 與“旅行者2號(hào)”的比較

“旅行者2號(hào)”最終在2007年距日84 AU處穿越終止激波[2]，之后一直在內(nèi)日球?qū)忧蕝^(qū)直至2018年底被確認(rèn)穿越日球?qū)禹斶M(jìn)入星際空間[23]。給出2011—2017年間的模擬與觀測(cè)的動(dòng)壓對(duì)比，在此期間“旅行者2號(hào)”從距日94.17 AU處運(yùn)動(dòng)至113.14 AU，如圖2所示。綠色橫線(xiàn)段表示發(fā)生在星際空間的等離子體振蕩事件，由“旅行者1號(hào)”攜帶的等離子體電場(chǎng)波動(dòng)儀器測(cè)量得到[24]。綠色的豎線(xiàn)對(duì)應(yīng)于“旅行者1號(hào)”磁強(qiáng)計(jì)測(cè)量得到的磁場(chǎng)跳變信號(hào)，對(duì)應(yīng)于星際空間傳播的激波[9]。

一般認(rèn)為圖2中的前3個(gè)激波信號(hào)與等離子體振蕩時(shí)間關(guān)聯(lián)得較好，后者被認(rèn)為來(lái)自星際激波前方的低能電子束流，類(lèi)似于II型的太陽(yáng)射電爆發(fā)信號(hào)[25]。第4個(gè)等離子體振蕩信號(hào)由于沒(méi)有對(duì)應(yīng)的磁場(chǎng)跳變，一般認(rèn)為它可能來(lái)源于星際其他地方的激波且未被“旅行者1號(hào)”的磁強(qiáng)計(jì)探測(cè)得到。6個(gè)字母A～F和對(duì)應(yīng)的垂直虛線(xiàn)表示“旅行者2號(hào)”觀測(cè)到的太陽(yáng)風(fēng)合并作用區(qū)動(dòng)壓的局地最高點(diǎn)，這些壓力脈沖被認(rèn)為與“旅行者1號(hào)”探測(cè)到的等離子體爆發(fā)事件信號(hào)相關(guān)，其中C、D和E被認(rèn)為對(duì)應(yīng)于第1、3和4個(gè)等離子體振蕩[15]。

圖2 2011—2017年采用3種不同內(nèi)邊界數(shù)據(jù)源的太陽(yáng)風(fēng)動(dòng)壓對(duì)比Fig.2 The comparision of solar wind dynamical pressure between the observation（red）and the simulations（blue）during the years 2011—2017

由于沒(méi)有終止激波，模擬得到的動(dòng)壓相對(duì)觀測(cè)值較大；出于方便考慮，將觀測(cè)的動(dòng)壓結(jié)果乘以1.5后與模擬結(jié)果做對(duì)比。圖2中可見(jiàn)，很難將模擬結(jié)果與“旅行者2號(hào)”的觀測(cè)做精確的逐一比較，一方面是模擬的有效區(qū)域在近赤道面，與“旅行者2號(hào)”的緯度存在較大差異，另外球?qū)ΨQ(chēng)流不存在與徑向方向相切的流動(dòng)，實(shí)際的太陽(yáng)風(fēng)傳播要復(fù)雜得多，此外終止激波的缺失也是一方面，以及其他的物理因素如湍流、拾起離子的影響等。

大體上，3個(gè)不同的內(nèi)邊界輸入源會(huì)導(dǎo)致不同的模擬結(jié)果，但趨勢(shì)大致相近，存在近似1年左右的長(zhǎng)周期變化[26]。觀測(cè)數(shù)據(jù)發(fā)現(xiàn)，激波穿越地球弓激波后速度為原先速度的0.7～1倍左右[27]。因此，設(shè)終止激波在84 AU左右，激波速度在～500 km/s，假設(shè)其在圖2中距離（94～113 AU）上穿越終止激波后以350 km/s運(yùn)動(dòng)10～30 AU左右，預(yù)計(jì)圖2中的模擬激波信號(hào)比實(shí)際傳播快0.04～0.12年左右，這種差異在圖中的刻度范圍內(nèi)并不顯著。從圖2中看，OMNI源的模擬結(jié)果對(duì)A、B、C和D動(dòng)壓脈沖具有一定的表現(xiàn)，但對(duì)后期兩個(gè)壓力脈沖E和F并不明顯；而STA源的模擬結(jié)果對(duì)B、D和F具有一定的表現(xiàn)，對(duì)A、C和E表現(xiàn)不明顯；對(duì)STB源的模擬結(jié)果來(lái)說(shuō)，前期在動(dòng)壓的變化趨勢(shì)上比較一致，但對(duì)后期的脈沖比如F表現(xiàn)不明顯，這種情況很可能是因?yàn)镾TB數(shù)據(jù)在自2014年10月后的缺失引起。

模擬結(jié)果表明，對(duì)一維模擬來(lái)說(shuō)，1 AU處不同經(jīng)度上的內(nèi)邊界輸入在外日球?qū)訁^(qū)雖然在大尺度結(jié)構(gòu)上具有一定的一致性，但細(xì)節(jié)方面仍具有差異，原因在于OMNI，STA和STB在1 AU軌道上的存在經(jīng)度差異，它們對(duì)太陽(yáng)風(fēng)事件的測(cè)量并不總是一致，比如存在有些太陽(yáng)風(fēng)事件只被某一經(jīng)度上的飛船探測(cè)到的情況[28]。因此在與“旅行者2號(hào)”的觀測(cè)數(shù)據(jù)對(duì)比時(shí)，需要綜合考慮3個(gè)不同內(nèi)邊界輸入的情況。以往人們?cè)谀M外日球?qū)犹?yáng)風(fēng)傳播時(shí)側(cè)重于使用OMNI數(shù)據(jù)[13]，該模擬算例表明在更高維度上如何處理好內(nèi)邊界問(wèn)題是一個(gè)很具有挑戰(zhàn)性的難題。

2.2 與“旅行者1號(hào)”的比較

“旅行者1號(hào)”在2013年8月于距日121.6 AU處穿越了日球?qū)禹擺29]，自那以后人們開(kāi)始能夠持續(xù)測(cè)量星際空間的等離子體波動(dòng)[10]，銀河宇宙線(xiàn)[30]和磁場(chǎng)信號(hào)[9]。然而，由于“旅行者1號(hào)”的等離子體儀器自1980年開(kāi)始就已經(jīng)失效，無(wú)法將模擬得到的等離子體參數(shù)與觀測(cè)數(shù)據(jù)直接比較。

圖3 2012—2016年用3種不同內(nèi)邊界數(shù)據(jù)源，模擬得到的“旅行者1號(hào)”可能探測(cè)的動(dòng)壓數(shù)據(jù)Fig.3 The simulated dynamical pressures that could be encountered by Voyager 1 during the years 2012—2016

圖3分別給出了采用OMNI，STA和STB源作為內(nèi)邊界輸入模擬得到“旅行者1號(hào)”自2012—2016年間測(cè)量得到的等離子體動(dòng)壓變化。黑色豎虛線(xiàn)表示“旅行者1號(hào)”在2 012.65時(shí)刻觀測(cè)到的穿越日球?qū)禹數(shù)奈恢?，另?個(gè)豎線(xiàn)表征時(shí)刻2 012.904、2 013.356和2 014.644，分別對(duì)應(yīng)于“旅行者1號(hào)”通過(guò)磁場(chǎng)強(qiáng)度跳變觀測(cè)得到的前向激波FS1，后向激波RS和前向激波FS2[9]?？梢钥闯?，模擬得到3個(gè)激波附近均具有持續(xù)時(shí)間較長(zhǎng)的壓力脈沖結(jié)構(gòu)甚至激波組合，時(shí)間范圍大致包括2 012.4—2 013.0、2 013.2—2 013.8和2 014.6—2 015.2。

由于模擬并未包含終止激波和日球?qū)禹?，?shí)際上這些壓力脈沖在與終止激波和日球?qū)禹斪饔煤竽軌蜻M(jìn)一步形成激波，進(jìn)而在星際空間被探測(cè)到。以O(shè)MNI數(shù)據(jù)的第2個(gè)動(dòng)壓脈沖結(jié)構(gòu)中的2 013.36處的激波為例，注意到在模擬后的第1 086 d到達(dá)85 AU，而在第1 212 d到達(dá)120 AU左右（本文未列出），可以估算激波速度為482 km/s。如果激波在85 AU左右穿越終止激波后速度下降為70%，則在120 AU處激波會(huì)遲54 d左右到達(dá)。而在星際空間里，磁場(chǎng)觀測(cè)顯示快磁聲波速大概為40 km/s[31]，因此激波穿越日球?qū)禹敽笏俣葧?huì)降至原來(lái)的1/10左右，這樣越往后模擬到的動(dòng)壓脈沖實(shí)際對(duì)應(yīng)的激波信號(hào)延遲越大。Richardson等[15]提出經(jīng)過(guò)修正后，“旅行者2號(hào)”觀測(cè)到的動(dòng)壓脈沖C、D、E分別對(duì)應(yīng)第1、3和4個(gè)等離子體振蕩信號(hào)，而第2個(gè)振蕩信號(hào)來(lái)源仍然未知。比較圖2和圖3可見(jiàn)，圖2中的C，D和E基本對(duì)應(yīng)圖3中的三個(gè)動(dòng)壓脈沖組合，即模擬得到的3個(gè)動(dòng)壓脈沖也同樣可能對(duì)應(yīng)于實(shí)際觀測(cè)到的第1、3和4個(gè)等離子體振蕩信號(hào)。此外，從模擬結(jié)果看，仍然無(wú)法對(duì)第2個(gè)振蕩信號(hào)的來(lái)源做出判斷。

以往的三維數(shù)值模擬結(jié)果認(rèn)為第1前向激波FS1來(lái)自于太陽(yáng)風(fēng)的共轉(zhuǎn)作用區(qū)在外日球?qū)訁^(qū)的合并發(fā)展，第2個(gè)前向激波是由共轉(zhuǎn)作用區(qū)與行星際日冕物質(zhì)拋射共同作用的結(jié)果[13]?？赡艿暮笙蚣げ≧S最初被認(rèn)為與2012年3月爆發(fā)的大規(guī)模日冕物質(zhì)拋射事件有關(guān)[11]。

由于模型的限制，這里沒(méi)有對(duì)模擬得到的3個(gè)動(dòng)壓脈沖結(jié)構(gòu)的驅(qū)動(dòng)源做分析，后續(xù)將采用更為全面的日球?qū)尤驍?shù)值模擬開(kāi)展相關(guān)研究。

3 結(jié) 論

本文通過(guò)簡(jiǎn)化的一維太陽(yáng)風(fēng)?星際中性原子相互作用的磁流體力學(xué)數(shù)值模型，對(duì)太陽(yáng)風(fēng)在外日球?qū)訁^(qū)的演化進(jìn)行了數(shù)值模擬，其中內(nèi)邊界設(shè)在距日1 AU處，外邊界在200 AU的星際空間。采用1 AU軌道上不同經(jīng)度的OMNI，STA和STB這3個(gè)不同觀測(cè)數(shù)據(jù)作為輸入條件，時(shí)間范圍從2010—2017年。模擬結(jié)果表明，采用不同內(nèi)邊界輸入后，雖然整體上看太陽(yáng)風(fēng)在外日球?qū)訁^(qū)的大尺度結(jié)構(gòu)比較類(lèi)似，但細(xì)節(jié)表現(xiàn)會(huì)有不同，因此在做數(shù)值模擬時(shí)需要做綜合考慮。模擬了“旅行者1號(hào)”在星際空間可能探測(cè)到的大尺度動(dòng)壓脈沖結(jié)構(gòu)，這些結(jié)構(gòu)被認(rèn)為能夠與日球?qū)禹斚嗷プ饔貌⑦M(jìn)一步產(chǎn)生激波，從而觸發(fā)“旅行者1號(hào)”觀測(cè)到的等離子電場(chǎng)振蕩信號(hào)，因此可以認(rèn)為這類(lèi)結(jié)構(gòu)也與觀測(cè)到的星際激波信號(hào)相關(guān)。