丁嘉威， 王正濤,2*， 詹文臻， 劉美琴

1 武漢大學(xué)測繪學(xué)院， 武漢 430079 2 武漢大學(xué)地球空間環(huán)境與大地測量教育部重點(diǎn)實(shí)驗(yàn)室， 武漢 430079

0 引言

行星的磁流體內(nèi)核的流體動(dòng)力學(xué)模擬主要是通過求解包括科里奧利力在內(nèi)的Navier-Stokes方程，并結(jié)合Boussinesq近似下的磁流體動(dòng)力學(xué)(MHD)感應(yīng)方程、溫度(或熵)方程與化學(xué)成分(本文未考慮)方程來實(shí)現(xiàn)的，得益于數(shù)值算法的不斷改進(jìn)以及高性能計(jì)算機(jī)的發(fā)展，使得當(dāng)前發(fā)電機(jī)模型的數(shù)值模擬計(jì)算速度越來越快，描述各物理量的參數(shù)越來越多且日益接近真實(shí)值.

1995年，超級(jí)計(jì)算機(jī)Cray C90上成功模擬了一個(gè)包含有限導(dǎo)電內(nèi)核的三維自恰地球發(fā)電機(jī)模型(Glatzmaier and Roberts, 1995a)，同年又模擬了長達(dá)4萬年的地球磁場，并成功實(shí)現(xiàn)了磁場的翻轉(zhuǎn)(Glatzmaier and Roberts, 1995b).此后相關(guān)研究者們?cè)谝陨匣A(chǔ)上對(duì)發(fā)電機(jī)模型中的物理參數(shù)進(jìn)行修改或增加、或使用更為合理的物理參數(shù)及數(shù)值方法來得到新的發(fā)電機(jī)模型并進(jìn)行系統(tǒng)性的研究(Christensen et al., 1998; Christensen et al., 1999; Olson et al., 1999; Christensen et al., 2001; Kutzner and Christensen, 2002; Takahashi et al., 2007; Takahashi et al., 2008; Wicht, 2014; Matsui et al., 2016).地球內(nèi)核的超速旋轉(zhuǎn)現(xiàn)象的發(fā)現(xiàn)(Song and Richards, 1996)，利用地震得到的地球內(nèi)核表面的鑲嵌結(jié)構(gòu)(Krasnoshchekov et al., 2005)，地核條件下鐵、固態(tài)鐵和鐵硅混合物的導(dǎo)熱性和導(dǎo)電性(Pozzo et al., 2012, 2014)，利用高壓試驗(yàn)和三維x射線顯微層析成像技術(shù)推斷的核-幔邊界處的溫度數(shù)據(jù)(Nomura et al., 2014)，通過數(shù)值模擬分析得到的深內(nèi)部地球結(jié)構(gòu)對(duì)內(nèi)核平動(dòng)振蕩本征周期的影響(江穎等, 2014)，關(guān)于核幔耦合對(duì)地球自由核章動(dòng)的激發(fā)影響(崔小明等, 2018)等各類更精細(xì)的地球內(nèi)部數(shù)據(jù)或相關(guān)研究都為更真實(shí)的地球(磁流體)發(fā)電機(jī)模型提供了很好的數(shù)據(jù)條件.

對(duì)于內(nèi)核在發(fā)電機(jī)中的作用，Hollerbach和Jones(1993)關(guān)于內(nèi)核對(duì)地磁場波動(dòng)和翻轉(zhuǎn)的作用進(jìn)行了研究；Glatzmaier 和Roberts(1995b)提出了包含有限導(dǎo)電內(nèi)核的三維、自恰地球發(fā)電機(jī)模型，并成功實(shí)現(xiàn)了主磁場的翻轉(zhuǎn)；Wicht(2002)研究分析了有限導(dǎo)電內(nèi)核對(duì)地磁場翻轉(zhuǎn)頻率的影響，并詳細(xì)描述了引入內(nèi)核后的發(fā)電機(jī)方程；Reshetnyak(2016)研究了半徑不斷增大的內(nèi)核在發(fā)電機(jī)中的作用.而關(guān)于內(nèi)核在外核驅(qū)動(dòng)下的高轉(zhuǎn)速(通常每時(shí)間尺度達(dá)幾百弧度)、大波動(dòng)的問題，可以通過引入引力耦合來穩(wěn)定內(nèi)核的轉(zhuǎn)速(Aubert and Dumberry, 2011).需要注意的是此處的高轉(zhuǎn)速、大波動(dòng)是以慣性系為參考進(jìn)行描述的，而并非以外核系統(tǒng)作為參考.

本文以Takahashi等(2007, 2008)所提出的兩種數(shù)值發(fā)電機(jī)模型為基準(zhǔn)，研究有限導(dǎo)電內(nèi)核在兩種旋轉(zhuǎn)方式下對(duì)發(fā)電機(jī)模型的極向和環(huán)向動(dòng)能、極向和環(huán)向磁能、核-幔邊界處的磁場強(qiáng)度、磁雷諾數(shù)、磁極翻轉(zhuǎn)頻率、軸向偶極子強(qiáng)度、赤道對(duì)稱性、緯向性和磁通量集中度的影響，其中后四者用于研究類地發(fā)電機(jī)的成立條件時(shí)所用的四種參數(shù)(Christensen et al., 2010).

1 基本原理

通常，對(duì)于磁流體發(fā)電機(jī)模型的求解會(huì)從描述流體運(yùn)動(dòng)的Navier-Stokes方程出發(fā)：

(1)

最終得到其無量綱形式為：

(2)

其中p″為修正后的壓力項(xiàng)，ez為旋轉(zhuǎn)軸方向的單位矢量.無量綱參數(shù)的具體表達(dá)式分別為：

埃克曼數(shù)

(3)

瑞利數(shù)

(4)

普朗特?cái)?shù)

(5)

磁普朗特?cái)?shù)

(6)

式中，下標(biāo)“o”表示外核邊界處的相應(yīng)物理值，“i”表示內(nèi)核邊界處的相應(yīng)物理值，d為外核厚度，ΔT為內(nèi)外核溫度差，κ為熱擴(kuò)散率.

最后，將方程(2)與以下方程聯(lián)立：

(7)

(8)

(9)

(10)

即得到求解發(fā)電機(jī)模型的無量綱方程組.

將有限導(dǎo)電的內(nèi)核引入發(fā)電機(jī)模型后，內(nèi)核和外核的發(fā)電機(jī)方程就必須同時(shí)進(jìn)行求解，它們由磁場和水平電場的連續(xù)性來關(guān)聯(lián)：

B=Bi,r=ri;

EH=EHi,r=ri.

(11)

若內(nèi)核在外核的驅(qū)動(dòng)下發(fā)生轉(zhuǎn)動(dòng)，則其驅(qū)動(dòng)力包含洛倫茲力與黏性力.轉(zhuǎn)速可由角動(dòng)量方程進(jìn)行描述，其具體表達(dá)式為：

(12)

其中I是沿Z軸旋轉(zhuǎn)的內(nèi)核轉(zhuǎn)動(dòng)慣量，ωi為內(nèi)核旋轉(zhuǎn)角速度，ΓL、Γv分別對(duì)應(yīng)于洛倫茲力及黏性力.內(nèi)核為導(dǎo)電的情況下，其對(duì)應(yīng)的發(fā)電機(jī)方程為：

(13)

在以上基礎(chǔ)上通過更改地球內(nèi)核的電導(dǎo)率、旋轉(zhuǎn)方式來研究其對(duì)主磁場帶來的影響.

2 模型參數(shù)

2.1 基準(zhǔn)模型

本文選擇使用了Takahashi等(2007)用于研究磁極轉(zhuǎn)換所使用的模型及Takahashi等(2008)所提出的數(shù)值發(fā)電機(jī)中埃克曼數(shù)為2×10-5的模型作為基準(zhǔn)模型，分別簡稱為FT07和FT08，模型中的球諧展開階次為16階，徑向的切比雪夫多項(xiàng)式的展開項(xiàng)數(shù)為32，時(shí)間尺度方案選擇了黏性擴(kuò)散時(shí)間，其具體轉(zhuǎn)換方法即前文對(duì)時(shí)間t的變量代換.

模型FT07所對(duì)應(yīng)的物理參數(shù)為：E=1×10-5，Ra=8×107，Pm=0.5，Pr=1；模型FT08所對(duì)應(yīng)的物理參數(shù)為：E=2×10-5，Ra=2×107，Pm=1.5，Pr=1.

對(duì)于內(nèi)核電導(dǎo)率，基準(zhǔn)模型將其設(shè)定為絕緣，且與外核之間不存在任何相關(guān)物理作用的.

本文默認(rèn)所有內(nèi)核的密度均等于外核密度.

所有模型所對(duì)應(yīng)的內(nèi)核的機(jī)械邊界條件選擇使用剛性、不可穿透的邊界，這與真實(shí)的內(nèi)核-外核邊界物理上最為接近.

內(nèi)核的熱邊界條件默認(rèn)均為固定熵條件，即熵si=sbottom.

2.2 由外核驅(qū)動(dòng)旋轉(zhuǎn)的有限導(dǎo)電內(nèi)核模型(DOC, Driven by Outer Core)

包含有限導(dǎo)電內(nèi)核的模型的無量綱參數(shù)、內(nèi)核密度、內(nèi)核的機(jī)械邊界條件以及熱邊界條件均與進(jìn)行比對(duì)的相應(yīng)基準(zhǔn)模型保持一致，且內(nèi)核的電導(dǎo)率均認(rèn)為與外核相同(目前內(nèi)核的增長通常被認(rèn)為是外核物質(zhì)在內(nèi)核邊界上凝固的結(jié)果，因此認(rèn)為內(nèi)核電導(dǎo)率和外核電導(dǎo)率相同在物理上解釋可能更為合理).

此模型下，內(nèi)核旋轉(zhuǎn)被認(rèn)為是洛倫茲力和黏滯力共同作用下的結(jié)果.

2.3 固定旋轉(zhuǎn)速度的有限導(dǎo)電內(nèi)核的模型I(FSR I, Fixed Super Rotation I)

該模型的轉(zhuǎn)速使用對(duì)應(yīng)的DOC模型在相應(yīng)時(shí)間尺度內(nèi)的最大轉(zhuǎn)速.分別對(duì)應(yīng)于351.06rad/時(shí)間尺度和383.96rad/時(shí)間尺度.

需要注意的是，此模型下的內(nèi)核不會(huì)與外核存在任何有關(guān)力的作用，即認(rèn)為內(nèi)核是自發(fā)旋轉(zhuǎn)、不受外核影響的.

2.4 固定旋轉(zhuǎn)速度的有限導(dǎo)電內(nèi)核的模型II(FSR II, Fixed Super Rotation II)

此模型與FSR I的區(qū)別僅在于內(nèi)核的轉(zhuǎn)速不同，所使用的轉(zhuǎn)速為所對(duì)應(yīng)的DOC模型在對(duì)應(yīng)時(shí)間尺度下的內(nèi)核平均轉(zhuǎn)速，分別為：32.90rad/時(shí)間尺度和-22.68rad/時(shí)間尺度.

2.5 固定旋轉(zhuǎn)速度的有限導(dǎo)電內(nèi)核的模型III(FSR III, Fixed Super Rotation III)

與FSR II類似，其轉(zhuǎn)速對(duì)應(yīng)為：-296.69rad/時(shí)間尺度和-363.98rad/時(shí)間尺度.

3 模型間的比較分析

3.1 核-幔邊界處的磁場強(qiáng)度

核-幔邊界處的磁場強(qiáng)度通常使用埃爾澤塞爾數(shù)進(jìn)行描述：

(14)

其大小等于無量綱磁場強(qiáng)度的平方，該數(shù)也可以用于描述洛倫茲力與科里奧利力的比值.式中下標(biāo)rms表均方根.

根據(jù)圖1可以看出模型FT07及其對(duì)應(yīng)的幾種比對(duì)模型的磁場強(qiáng)度在0.4個(gè)時(shí)間尺度后基本在0.3左右上下波動(dòng)，因此對(duì)于模型FT07，其所有的對(duì)比均使用0.4～3個(gè)時(shí)間尺度內(nèi)的數(shù)據(jù).同理根據(jù)圖2，F(xiàn)T08使用0.5～3個(gè)時(shí)間尺度內(nèi)的數(shù)據(jù).

圖1 模型FT07及其對(duì)應(yīng)的DOC、FSR I、FSR II、FSR III的核-幔邊界處的磁場強(qiáng)度.該模型所使用的積分時(shí)間步長的數(shù)量級(jí)達(dá)10-6，造成了原圖曲線存在大量重疊的情況，因此此處每隔5000個(gè)點(diǎn)取一個(gè)來繪制其趨勢，計(jì)算時(shí)仍采用全部數(shù)據(jù).后作類似處理的圖像會(huì)標(biāo)注相應(yīng)的間隔，如“(5000)”Fig.1 Magnetic field intensity at core-mantle boundary of benchmark dynamo FT07 and its corresponding dynamo DOC, FSR I, FSR II, and FSR III. Order of magnitude of integral time step used in this model is up to 10-6, which causes a large amount of overlap in the original graph curve. Therefore, here every 5000 points are taken to plot its trend, but all data are still used in the calculation. Images processed similarly are marked with corresponding intervals, for example,“(5000)”

圖2 模型FT08及其對(duì)應(yīng)的DOC、FSR I、FSR II、FSR III的核-幔邊界處的磁場強(qiáng)度(5000)Fig.2 Magnetic field intensity at core-mantle boundary of benchmark dynamo FT08 and its corresponding dynamo DOC, FSR I, FSR II, and FSR III (5000)

3.2 DOC模型轉(zhuǎn)速

內(nèi)核轉(zhuǎn)速可由角動(dòng)量方程(12)進(jìn)行描述.計(jì)算得到的轉(zhuǎn)速圖見圖3—4.

圖3 模型FT07所對(duì)應(yīng)的DOC模型在0.01～3個(gè)時(shí)間尺度內(nèi)的內(nèi)核轉(zhuǎn)速(5000)Fig.3 Inner-core rotation speed of the DOC model corresponding to model FT07 in 0.01～3 timescales (5000)

圖4 模型FT08所對(duì)應(yīng)的DOC模型在0.01～3個(gè)時(shí)間尺度內(nèi)的內(nèi)核轉(zhuǎn)速(5000)Fig.4 Inner-core rotation speed of the DOC model corresponding to model FT08 in 0.01～3 timescales (5000)

3.3 外核能量

通過比較外核的動(dòng)能、磁能以及相應(yīng)的磁場變化來研究不同驅(qū)動(dòng)方式下的內(nèi)核對(duì)數(shù)值發(fā)電機(jī)作用的影響.在螺線向量場在球坐標(biāo)系下的分析中通常會(huì)對(duì)其應(yīng)用極向-環(huán)向分解，是亥姆霍茲分解的一個(gè)受限制的形式，在本文所涉及的問題中，動(dòng)能和磁能具體的極向-環(huán)向分解后表達(dá)形式分別為：

(15)

(16)

式中，Wl m為動(dòng)力場極向位能，Zl m為動(dòng)力場環(huán)向位能，bl m為磁場極向位能，jl m為磁場環(huán)向位能.

3.4 磁雷諾數(shù)

磁雷諾數(shù)的定義為：

(17)

主要用于描述磁流體內(nèi)部速度u和磁場B的方均根值的比值.

3.5 磁極翻轉(zhuǎn)頻率

磁極翻轉(zhuǎn)頻率f可根據(jù)磁極緯度的變化次數(shù)統(tǒng)計(jì)得到.

3.6 軸向偶極子強(qiáng)度

本文中軸向偶極子強(qiáng)度定義為偶極子(axial dipole,AD)與非偶極子(non-axial dipole,NAD)之比：

(18)

3.7 赤道對(duì)稱性

對(duì)于赤道對(duì)稱性，通常使用奇偶比O/E進(jìn)行描述，即2到8階的調(diào)和級(jí)數(shù)在核-幔邊界處的階次和(n+m)為奇數(shù)的分量和與階次和為偶數(shù)的分量和之比.

3.8 緯向性

緯向性即緯向-非緯向比Z/NZ是核-幔邊界處所有緯向性分量和非緯向性分量之比.

3.9 磁通量集中度

核-幔邊界處的磁通量在一定程度上集中并形成了強(qiáng)輻向場斑塊Br，可用通量集中因子FCF來定量描述：

(19)

其中〈〉表示在整個(gè)球面上求平均.

3.10 比較分析

表1給出了包含導(dǎo)電內(nèi)核的各發(fā)電機(jī)模型的各物理量較基準(zhǔn)模型的變化量.可以看出，有限導(dǎo)電內(nèi)核的引入對(duì)模型的磁場強(qiáng)度只會(huì)產(chǎn)生最高不到4%的變化，極、環(huán)向動(dòng)能最高僅約2%，極、環(huán)向磁能最高也約為4%，磁雷諾數(shù)的變化量不到1%，磁通量強(qiáng)度最高約5%.但對(duì)于軸向偶極子強(qiáng)度和磁極翻轉(zhuǎn)頻率，各模型相較基準(zhǔn)模型的變化量則存在較大差異，特別的，DOC模型對(duì)于磁極翻轉(zhuǎn)頻率幾乎沒有影響；對(duì)于赤道對(duì)稱性，除以FT07為基準(zhǔn)的FSR I的變化量達(dá)到+19.03%之外，其他的變化量均小于6%；在緯向性的比較中，除了以FT07為基準(zhǔn)的FSR II和FSR III的變化量達(dá)17%以上，其他模型的變化量均在10%以下.

表1 各模型對(duì)應(yīng)的各物理量較基準(zhǔn)模型的變化量Table 1 Variation of each physical quantity corresponding to each model compared with the benchmark dynamo

4 結(jié)論

有限導(dǎo)電內(nèi)核的引入，對(duì)于模型的磁場強(qiáng)度、極向和環(huán)向動(dòng)能、極向和環(huán)向磁能、磁雷諾數(shù)和磁通量集中度的影響是極弱的，只會(huì)產(chǎn)生小于6%的變化量.這可能是由于內(nèi)核本身的能量相對(duì)于外核來說是極小的，因此引入有限導(dǎo)電的內(nèi)核后，其所帶來的影響也應(yīng)不會(huì)明顯.

圖5 模型FT07及其對(duì)應(yīng)的DOC、FSR I、FSR II、FSR III的磁極緯度Fig.5 Latitude of magnetic pole of benchmark dynamo FT07 and its corresponding dynamo DOC, FSR I, FSR II, and FSR III (5000)

對(duì)于軸向偶極子強(qiáng)度，由于本文所使用的模型是可以實(shí)現(xiàn)磁場翻轉(zhuǎn)的發(fā)電機(jī)模型，而這類模型的軸向偶極子強(qiáng)度通常都較弱，其本身基數(shù)較小，因此引入有限導(dǎo)電內(nèi)核后，導(dǎo)致其相對(duì)變化量較大，但其最大絕對(duì)變化量小于0.03，且結(jié)合圖5—8，可以明顯看出，有限導(dǎo)電內(nèi)核對(duì)軸向偶極子強(qiáng)度的影響更多的集中在磁極翻轉(zhuǎn)的時(shí)刻，這表明對(duì)于那些偶極性較弱的發(fā)電機(jī)模型，有限導(dǎo)電內(nèi)核的引入與否會(huì)直接影響該發(fā)電機(jī)的偶極子強(qiáng)度，并且對(duì)于可實(shí)現(xiàn)磁場翻轉(zhuǎn)的發(fā)電機(jī)模型來說，這一影響可能會(huì)更明顯.

而對(duì)于赤道對(duì)稱性和緯向性，除部分模型存在較大的變化外，其余模型的影響均小于10%，若發(fā)電機(jī)模型涉及以上兩個(gè)物理量，有限導(dǎo)電內(nèi)核(FSR)仍是值得考慮的因素.

對(duì)于翻轉(zhuǎn)頻率來說，DOC模型的翻轉(zhuǎn)次數(shù)與基準(zhǔn)模型幾乎完全一致，而FSR模型則對(duì)翻轉(zhuǎn)頻率造成了不同程度的較明顯的影響.這表明DOC模型的內(nèi)核旋轉(zhuǎn)雖然波動(dòng)大、轉(zhuǎn)速高，但可能由于其驅(qū)動(dòng)

圖6 模型FT08及其對(duì)應(yīng)的DOC、FSR I、FSR II、FSR III的磁極緯度Fig.6 Latitude of magnetic pole of benchmark dynamo FT08 and its corresponding dynamo DOC, FSR I, FSR II, and FSR III (5000)

圖7 模型FT07及其對(duì)應(yīng)的DOC、FSR I、FSR II、FSR III的軸向偶極子強(qiáng)度(5000)Fig.7 Relative axial dipole power of benchmark dynamo FT07 and its corresponding dynamo DOC, FSR I, FSR II, and FSR III (5000)

圖8 模型FT08及其對(duì)應(yīng)的DOC、FSR I、FSR II、FSR III的軸向偶極子強(qiáng)度(5000)Fig.8 Relative axial dipole power of benchmark dynamo FT08 and its corresponding dynamo DOC, FSR I, FSR II, and FSR III (5000)

方式及機(jī)理與真實(shí)情況相近且考慮到內(nèi)核本身的能量對(duì)外核產(chǎn)生的影響程度有限，因此DOC模型與基準(zhǔn)模型的翻轉(zhuǎn)頻率幾乎完全一致也在物理上更為合理；而對(duì)于FSR模型，其本身就是一種為了方便控制轉(zhuǎn)速而理想化后的模型，因此其內(nèi)核可能會(huì)對(duì)外核磁流體產(chǎn)生一定影響，而引起磁極翻轉(zhuǎn)頻率發(fā)生較明顯變化.

總體來說，本文中的DOC模型除偶極性以外，對(duì)于其他物理量所能帶來的變化都是可忽略的，且考慮到實(shí)際內(nèi)核較外核的能量比來說，以上結(jié)果是較合理的.

致謝MagIC開發(fā)團(tuán)隊(duì)提供的球殼中的流體動(dòng)力學(xué)數(shù)值代碼MagIC、南方科技大學(xué)林玉峰副教授提供的理論指導(dǎo)及武漢大學(xué)超算中心提供的算力支持.