韓宇峰馬紹賢蘇彩虹

(天津大學(xué)機(jī)械學(xué)院 高速空氣動(dòng)力學(xué)研究室,天津 300072)

0 引 言

高超聲速三維邊界層的轉(zhuǎn)捩預(yù)測(cè)是飛行器設(shè)計(jì)中必須考慮的重要問題,特別是對(duì)長(zhǎng)時(shí)間在大氣中飛行的高速飛行器。由于湍流狀態(tài)下壁面摩阻和熱流比層流下大得多,只有準(zhǔn)確地預(yù)測(cè)出轉(zhuǎn)捩位置,才能準(zhǔn)確計(jì)算阻力和表面熱流,適當(dāng)設(shè)計(jì)熱防護(hù)措施。實(shí)際中飛行器飛行時(shí)常常有后掠角,這時(shí)在邊界層內(nèi)與外流流線垂直的方向上由于壓力梯度會(huì)產(chǎn)生二次流動(dòng),稱為橫流。橫流失穩(wěn)是觸發(fā)三維邊界層轉(zhuǎn)捩至湍流的主要原因。然而,目前對(duì)橫流失穩(wěn)觸發(fā)轉(zhuǎn)捩發(fā)生的機(jī)理,特別是對(duì)高超聲速三維邊界層[1-4],還并不清楚。

關(guān)于橫流不穩(wěn)定性的研究,早期主要針對(duì)低速(不可壓或亞聲速、個(gè)別低超聲速)后掠翼和后掠平板的邊界層。國(guó)際上最具代表性的是德宇航Bippes及其合作者、以及美國(guó)Saric及其合作者的工作,比較全面的綜述可見文獻(xiàn)[3-4]。國(guó)內(nèi)也開展了不少關(guān)于橫流感受性[5]、橫流穩(wěn)定性[6-8]以及橫流轉(zhuǎn)捩預(yù)測(cè)模型的研究[9-10]。目前,已經(jīng)有了比較一致的認(rèn)識(shí),即有橫流的邊界層中可以存在兩種由橫流不穩(wěn)定性導(dǎo)致的流動(dòng)——定常橫流渦和低頻的橫流行進(jìn)波。定常渦對(duì)壁面粗糙度非常敏感,很容易被激發(fā)。激發(fā)出的定常渦經(jīng)歷一定程度的線性增長(zhǎng)后,在下游會(huì)由于非線性作用幅值出現(xiàn)飽和。飽和的定常渦導(dǎo)致平均流發(fā)生修正,使得在修正的平均流上高頻擾動(dòng)再次失穩(wěn),即出現(xiàn)二次失穩(wěn)現(xiàn)象,促使轉(zhuǎn)捩發(fā)生[11]。橫流轉(zhuǎn)捩中定常渦和行進(jìn)波的非線性飽和過程很長(zhǎng),在此過程中,擾動(dòng)幅值變化非常緩慢,轉(zhuǎn)捩發(fā)生的位置與開始發(fā)生非線性飽和的位置有相當(dāng)長(zhǎng)的距離。因此,若仍采用目前普遍用于Mack模態(tài)轉(zhuǎn)捩的預(yù)測(cè)方法,即以首次不穩(wěn)定波(橫流定常渦或行進(jìn)波)的幅值為基礎(chǔ)作為轉(zhuǎn)捩判據(jù)是不可行的[11-14]。若以此為判據(jù),任何在輸入擾動(dòng)和計(jì)算擾動(dòng)演化過程中引入的很小的誤差,也會(huì)給轉(zhuǎn)捩位置的預(yù)測(cè)帶來很大的影響。因此,橫流失穩(wěn)導(dǎo)致的轉(zhuǎn)捩預(yù)測(cè)相比Mack模態(tài)失穩(wěn)而言,更依賴于對(duì)轉(zhuǎn)捩機(jī)理的認(rèn)識(shí)和理解。

本文針對(duì)高超聲速后掠鈍板的三維邊界層,采用直接數(shù)值模擬方法(DNS)和全局不穩(wěn)定性(Biglobal)的分析方法,研究了橫流定常渦模態(tài)的首次失穩(wěn)和二次失穩(wěn)；采用DNS計(jì)算了高超聲速三維邊界層橫流失穩(wěn)轉(zhuǎn)捩至湍流的過程,分析了邊界層轉(zhuǎn)捩過程的機(jī)理。

1 基本流場(chǎng)

1.1 計(jì)算工況

本文選取半無限長(zhǎng)后掠鈍板為研究對(duì)象,其前緣為一半徑為3.5 cm的圓柱形,如圖1所示。x、y、z分別為流向、法向、展向。飛行馬赫數(shù)為6,來流迎角為0°,后掠角為45°。飛行高度為30 km,來流溫度為226.5 K,以頭部半徑為特征長(zhǎng)度的雷諾數(shù)為7.91×104。

圖1 后掠鈍板模型示意圖Fig.1 Schematic diagram of the swept blunt plate

1.2 計(jì)算方法和邊界條件

直角坐標(biāo)系下三維N-S方程的守恒型控制方程為:

其中:U為守恒型變量分別表示密度、流向速度、法向速度、展向速度和溫度,c v為定容比熱；E、F、G分別為流向、法向和展向的無黏通量；E v、F v、G v為相應(yīng)的黏性通量。N-S方程中具體系數(shù)矩陣見文獻(xiàn)[20]。

黏性系數(shù)采用Sutherland公式計(jì)算,具體為μ＝下標(biāo)∞表示來流量。熱傳導(dǎo)系數(shù)κ由給出,Pr＝0.72,c p∞為來流定壓比熱。

“中國(guó)和印度的經(jīng)驗(yàn)表明，強(qiáng)大國(guó)家和強(qiáng)大社會(huì)同時(shí)出現(xiàn)，隨著時(shí)間的流逝，而相互平衡，相互抵消，這樣才會(huì)有較好形式的自由”。③馬克思·韋伯認(rèn)為，中國(guó)在秦朝就具有現(xiàn)代國(guó)家治理的要素，但是沒有形成強(qiáng)大社會(huì)相平衡，法治要素缺乏，所以難以形成現(xiàn)代民主化的國(guó)家構(gòu)建；相反，印度則因宗教、種姓制度形成了強(qiáng)大的社會(huì)要素基礎(chǔ)，但難以形成強(qiáng)大的統(tǒng)一國(guó)家，也影響了現(xiàn)代國(guó)家的構(gòu)建。英國(guó)基于個(gè)人主義、權(quán)利制衡基礎(chǔ)上的政治治理結(jié)構(gòu)，是特定條件下形成的特定結(jié)果。

由于模型是二維的,展向無限長(zhǎng),因此?/?z＝0。但由于流動(dòng)存在后掠角,因此是三維的,即速度分量包含u、v和w。

采用基于有限差分的直接數(shù)值模擬程序計(jì)算基本流場(chǎng)。對(duì)流項(xiàng)離散采用五階WENO格式,在邊界處需進(jìn)行降階處理[22]。對(duì)于邊界點(diǎn),采用三點(diǎn)二階的偏心差分格式；對(duì)于次邊界點(diǎn),采用四點(diǎn)三階的偏心差分格式。黏性項(xiàng)離散采用六階中心差分格式,時(shí)間推進(jìn)采用四步四階Runge-Kutta法。壁面采用無滑移和絕熱邊界條件,遠(yuǎn)場(chǎng)條件給定自由來流。頭部對(duì)稱面處采用對(duì)稱邊界條件,出口采用線性外推邊界條件。本文采用的基本流程序代碼已成功地用于后掠鈍板的基本流計(jì)算中,可見文獻(xiàn)[21]。

由于橫流渦擾動(dòng)和二次失穩(wěn)擾動(dòng)是三維的,因此,對(duì)于定常橫流渦的演化,以及非定常二次失穩(wěn)模態(tài)的演化,均采用三維N-S方程。

擾動(dòng)演化的計(jì)算中,入口給定具體的擾動(dòng)波形式,出口采用嵌邊區(qū)[23],壁面采用無滑移邊界條件,外邊界采用遠(yuǎn)場(chǎng)邊界條件,展向采用周期邊界條件。數(shù)值離散格式采用與基本流計(jì)算中相同的格式。

1.3 基本流場(chǎng)

計(jì)算域如圖2所示。流向計(jì)算域?yàn)?00個(gè)頭部半徑(3.5m),法向計(jì)算域包含激波。流向和法向網(wǎng)格數(shù)分別為1001和301,CFL數(shù)取為0.4。

圖2 基本流計(jì)算域示意圖Fig.2 Schematic diagram of the computing domain

求解二維的N-S方程直至所有網(wǎng)格點(diǎn)上的變量達(dá)到定常狀態(tài)。將速度在邊界層外緣的勢(shì)流方向和垂直于勢(shì)流的方向進(jìn)行分解,得到勢(shì)流方向速度分量u s和橫流速度分量u c的分布,如圖3所示?？梢钥吹?流動(dòng)中橫流的量級(jí)為10-2。從頭部到下游,橫流強(qiáng)度逐漸減弱。

2 首次失穩(wěn)——橫流定常渦的演化

首先對(duì)基本流進(jìn)行線性穩(wěn)定性(LST)分析。在平行流假設(shè)下,小擾動(dòng)可以寫成如下形式:

其中(y)為擾動(dòng)的形函數(shù),擾動(dòng)在流向和展向和時(shí)間方向具有波的形式?？紤]空間模式,即擾動(dòng)沿空間增長(zhǎng),那么頻率ω為實(shí)數(shù),α和β為復(fù)數(shù),其虛部的相反數(shù)分別表示擾動(dòng)波在流向和展向的增長(zhǎng)率。將上式帶入到線性化的N-S方程的擾動(dòng)形式中,并忽略基本流量的流向?qū)?shù)及法向速度,得到準(zhǔn)平行流假設(shè)下的線性穩(wěn)定性方程:

在壁面和遠(yuǎn)場(chǎng)擾動(dòng)為零。這樣方程結(jié)合齊次邊界條件,構(gòu)成特征值問題。通過求解可以得到邊界層擾動(dòng)波的色散關(guān)系和擾動(dòng)的特征函數(shù)。

圖3 勢(shì)流方向速度(a)和橫流方向速度(b)Fig.3 Potential flow(a)and crossflow velocities(b)

橫流有兩種失穩(wěn)模態(tài),一種是定常的橫流渦,另一種是低頻的橫流行進(jìn)波。在高空低背景擾動(dòng)情況下,定常橫流渦被認(rèn)為是三維邊界層中占據(jù)主導(dǎo)作用的失穩(wěn)模態(tài)[3]。因此,在首次失穩(wěn)中,我們暫時(shí)只關(guān)注定常橫流渦。本文選取β0＝3這一定常渦擾動(dòng),來計(jì)算其在流場(chǎng)中的演化情況。

由于基本流中引入的是三維擾動(dòng),因此首先將基本流在展向延拓成一個(gè)波長(zhǎng),即L z＝2π/β0,然后在計(jì)算域的入口引入擾動(dòng)來研究首次失穩(wěn)擾動(dòng)演化,擾動(dòng)形式為:

其中,A0代表入口擾動(dòng)幅值,取為10－3(y)代表線性穩(wěn)定性分析得到的β0擾動(dòng)的特征函數(shù)。

用DNS計(jì)算這一擾動(dòng)在流場(chǎng)中的演化,直至流場(chǎng)定常。圖4給出了定常橫流渦的流向速度云圖?？梢钥闯鲈谏嫌?流向速度在展向分布比較均勻,隨著擾動(dòng)向下游演化,壁面附近的流體逐漸被卷起,形成渦結(jié)構(gòu)。

將DNS得到的擾動(dòng)場(chǎng)沿展向做傅里葉變換,可以得到不同展向波數(shù)的各階譜的速度擾動(dòng)幅值沿流向的變化,如圖5所示。圖中(0,1)表示加入的定?；静?β0＝3),(0,2)和(0,3)分別表示2倍展向波數(shù)和3倍展向波數(shù)的波?？梢钥闯?在上游基本波的幅值快速增長(zhǎng),在x＝20處由于幅值較大開始非線性作用產(chǎn)生基本流修正(0,0)和二次諧波(0,2)。當(dāng)基本流修正增長(zhǎng)到幅值約為0.1左右時(shí),(0,1)波幅值達(dá)到飽和,隨后基本流修正也開始出現(xiàn)飽和。再向下游演化,擾動(dòng)幅值不再增大,擾動(dòng)維持在接近飽和的狀態(tài)。

圖4 定常橫流渦流向速度云圖Fig.4 Contours of the streamwise velocity at different locations

圖5 擾動(dòng)各階譜的流向速度幅值Fig.5 The modal amplitudes of the streamwise fluctuation velocity

由于定常渦的非線性作用會(huì)產(chǎn)生基本流修正,修正后基本流的失穩(wěn)特征將會(huì)出現(xiàn)變化。因此我們將圖5中得到的基本流修正項(xiàng)(0,0)加入到基本流中,再次分析穩(wěn)定性特征。圖6給出了新的基本流在不同流向位置的失穩(wěn)區(qū)內(nèi)增長(zhǎng)率等值線圖。同時(shí),作為對(duì)比,也給出了原始基本流的穩(wěn)定性分析結(jié)果?？梢钥闯?在x＝20的地方,新的基本流增長(zhǎng)率分布與原始基本流一致,這是由于基本流修正的幅值還比較小,對(duì)穩(wěn)定性影響不大。在下游x＝30的位置,相比原始基本流,增長(zhǎng)率以及失穩(wěn)區(qū)域明顯減??；繼續(xù)向下游,在x＝40處新的基本流幾乎不再失穩(wěn)。說明由于基本流修正的作用,流場(chǎng)不穩(wěn)定區(qū)域減小,這意味著忽略擾動(dòng)展向的變化,只考慮(0,0)基本流修正,只會(huì)得到更穩(wěn)定的基本流。但是轉(zhuǎn)捩的發(fā)生依賴于流場(chǎng)在更大范圍內(nèi)的失穩(wěn),說明具有展向尺度的二次失穩(wěn)擾動(dòng)可能是促進(jìn)轉(zhuǎn)捩的關(guān)鍵因素。

圖6 線性穩(wěn)定性分析得到的增長(zhǎng)率等值線圖Fig.6 Growth rate contours obtained by linear stability analysis(a)for the mean flow profile after primary instability and(b)the base flow profile

3 二次失穩(wěn)及轉(zhuǎn)捩至湍流的過程

首先對(duì)首次失穩(wěn)的流場(chǎng)進(jìn)行二維全局穩(wěn)定性分析。全局穩(wěn)定性分析方法可以研究基本流在兩個(gè)方向是快變而在另一個(gè)方向是慢變的三維流場(chǎng)的穩(wěn)定性問題。擾動(dòng)可以寫成如下形式:

其中(y,z)為擾動(dòng)的形函數(shù),擾動(dòng)在流向和時(shí)間方向具有波的形式。與LST類似,經(jīng)過推導(dǎo)可以得到關(guān)于(y,z)的二維特征值問題。具體公式見文獻(xiàn)[21]。

對(duì)首次失穩(wěn)后流場(chǎng)的某一個(gè)流向站位的剖面作全局穩(wěn)定性分析,可以找到不穩(wěn)定的模態(tài),稱為二次失穩(wěn)模態(tài)。

在x＝36處,即首次失穩(wěn)波接近飽和(幅值約為0.2)位置處進(jìn)行全局不穩(wěn)定分析,可以找到兩個(gè)二次模態(tài),分別對(duì)應(yīng)于z模態(tài)和y模態(tài)[13]。圖7(a)(b)分別給出了z模態(tài)和y模態(tài)的特征函數(shù)對(duì)于z模態(tài),其最大值的位置在法向y≈1.5,即二次失穩(wěn)擾動(dòng)發(fā)生在飽和渦結(jié)構(gòu)的側(cè)部,這一類型的模態(tài)是由流場(chǎng)中展向速度剪切導(dǎo)致的。而對(duì)于y模態(tài),其最大值的位置在法向y≈1.6,對(duì)應(yīng)飽和渦結(jié)構(gòu)的頂部位置,這一類型的模態(tài)是由流場(chǎng)中法向速度剪切導(dǎo)致的。圖8給出了x＝36處,時(shí)間模式的全局穩(wěn)定性分析無量綱頻率為2的y模態(tài)(b)的特征函數(shù)

圖7 x＝36,無量綱頻率為5的z模態(tài)(a)和

圖8 x＝36,二次失穩(wěn)模態(tài)的增長(zhǎng)率Fig.8 Growth rates of the secondary instability mode at x＝36

以x＝36為計(jì)算域入口,引入全局穩(wěn)定性分析得到的兩個(gè)不同頻率(ω＝5和ω＝4)的二次失穩(wěn)模態(tài)。為了使轉(zhuǎn)捩盡快發(fā)生,它們的初始幅值都選為A0＝0.01。使用DNS計(jì)算二次失穩(wěn)擾動(dòng)演化,其中流向計(jì)算域?yàn)閤＝36～96,展向計(jì)算域仍為一個(gè)基本波的波長(zhǎng)。流向、法向和展向計(jì)算域內(nèi)的網(wǎng)格點(diǎn)數(shù)分別為2401、151和61。

圖9給出了計(jì)算得到的流場(chǎng)瞬時(shí)速度等值面?？梢钥吹皆谏嫌慰拷肟诘奈恢?二次失穩(wěn)擾動(dòng)幅值很小,流場(chǎng)的渦結(jié)構(gòu)主要呈現(xiàn)首次失穩(wěn)后接近飽和的定常渦結(jié)構(gòu)。隨著擾動(dòng)向下游演化,高頻的二次失穩(wěn)擾動(dòng)增長(zhǎng)起來,飽和渦結(jié)構(gòu)開始扭曲,繼續(xù)向下游演化,規(guī)則的渦結(jié)構(gòu)破碎,轉(zhuǎn)捩完成,最終形成湍流。

圖9 瞬時(shí)流向速度等值面(0.55)圖Fig.9 Contours of the streamwise velocity at the value of 0.55

圖10給出了不同流向位置的二次失穩(wěn)擾動(dòng)速度云圖以及首次失穩(wěn)橫流渦速度等值線圖。從圖中可以看出在x＜51的范圍內(nèi),二次失穩(wěn)擾動(dòng)的形狀與首次失穩(wěn)橫流渦的速度等值線形狀相近,說明首次失穩(wěn)流場(chǎng)的剪切導(dǎo)致流動(dòng)出現(xiàn)二次失穩(wěn)。而在x＞51的下游,明顯觀察到二次失穩(wěn)擾動(dòng)形狀逐漸變得不規(guī)則,說明流場(chǎng)中不同頻率(展向波數(shù))的高次諧波增長(zhǎng)起來,流場(chǎng)開始出現(xiàn)轉(zhuǎn)捩。

為了進(jìn)一步分析流場(chǎng)的轉(zhuǎn)捩過程,首先分析二次失穩(wěn)擾動(dòng)各階譜的演化特征。圖11給出了傅里葉分析得到的不同頻率的擾動(dòng)幅值沿流向的變化,圖中ω＝4和ω＝5的波表示入口加入的兩個(gè)基本波,無量綱頻率ω＝0、1、2、3的波表示由基本波相互非線性作用激發(fā)的擾動(dòng),更高頻率的擾動(dòng)在圖中用灰線表示。可以看到,在x＜44的范圍內(nèi),入口加入的基本波快速增長(zhǎng),同時(shí)由于它們的幅值較大,非線性作用得到的ω＝0和ω＝1的擾動(dòng)被激發(fā)并且快速增長(zhǎng)。繼續(xù)向下游演化,其它頻率的擾動(dòng)也被激發(fā),并且它們都會(huì)經(jīng)歷一個(gè)快速增長(zhǎng)階段,最終在x＝55附近各階譜的幅值趨于穩(wěn)定,其中頻率為0的譜幅值最大,并且頻率越大,幅值越小。

圖10 不同流向位置的二次失穩(wěn)擾動(dòng)速度云圖以及首次失穩(wěn)橫流渦速度等值線圖Fig.10 Contours of the disturbance velocities for the secondary instability mode and the contour lines of the streamwise velocities for the cross flow vortex

圖11 傅里葉變換得到的各階譜的幅值Fig.11 Fourier modal of amplitude evolution

圖12給出了壁面摩擦系數(shù)(C f)沿流向的變化,其中紫線是層流解的結(jié)果,紅線是采用SST模型得到的湍流剖面的結(jié)果,而黃線是首次失穩(wěn)后橫流渦的壁面摩擦系數(shù),黑色實(shí)線是二次失穩(wěn)后的壁面摩擦系數(shù)?？梢钥吹绞状问Х€(wěn)后壁面摩擦系數(shù)有一定的抬升,但是抬升量不大,并沒有達(dá)到湍流的壁面摩擦系數(shù)值。說明首次失穩(wěn)幅值增大到一定程度對(duì)平均流產(chǎn)生修正,導(dǎo)致壁面摩阻上升,但是由于首次失穩(wěn)波已經(jīng)飽和,其本身無法再使得壁面剪切增大,使得C f進(jìn)一步抬升,因此轉(zhuǎn)捩并未發(fā)生。而加入二次失穩(wěn)擾動(dòng)之后,壁面摩擦系數(shù)在40＜x＜60的范圍內(nèi)快速抬升,并達(dá)到湍流剖面的摩擦系數(shù),說明在x＝40下游,流場(chǎng)開始轉(zhuǎn)捩,并逐步演化為湍流。這是由于引入高頻的二次失穩(wěn)波后,在非線性作用下,各階擾動(dòng)經(jīng)歷快速的增長(zhǎng)(見圖11),不斷從平均流中吸取能量,從而進(jìn)一步修正平均流使得C f更加快速抬升。這樣的結(jié)果,導(dǎo)致壁面摩擦系數(shù)超過了相應(yīng)的湍流的估計(jì)值,因而產(chǎn)生了“過沖”[24]的現(xiàn)象。當(dāng)各階擾動(dòng)的幅值達(dá)到一定程度,不再繼續(xù)增長(zhǎng),“過沖”現(xiàn)象逐漸減弱或消失,C f曲線會(huì)逐漸趨向于充分發(fā)展湍流的估計(jì)值。同時(shí)在C f曲線快速抬升的階段,低頻擾動(dòng)快速增長(zhǎng)并且幅值最大,說明低頻擾動(dòng)的增長(zhǎng)主導(dǎo)平均流修正及C f曲線的抬升。

圖12 壁面摩擦系數(shù)沿流向的變化Fig.12 Variation of the wall friction coefficient

將展向計(jì)算域拓展為兩個(gè)基本波波長(zhǎng)(2L z),其他條件保持不變,計(jì)算了二次失穩(wěn)擾動(dòng)的演化,圖12中黑色虛線給出了壁面摩擦系數(shù)(C f)沿流向的變化,與展向計(jì)算域?yàn)橐粋€(gè)基本波波長(zhǎng)(黑色實(shí)線)結(jié)果相比,二者幾乎一致,特別是在C f曲線抬升階段,二者完全重合,說明展向計(jì)算域?qū)M流渦轉(zhuǎn)捩過程影響不大。

圖13 平均速度剖面的分布Fig.13 Distribution of the mean velocity

圖13給出了經(jīng)過Van Driest變換后的平均速度剖面的分布,圖中虛線分別表示壁面律(u＋＝y(tǒng)＋)和對(duì)數(shù)律(u＋＝2.5lny＋＋4)?？梢钥闯鲈?1＜x＜91的范圍內(nèi),平均速度剖面與對(duì)數(shù)率和壁面率符合得很好,而且隨著向下游演化,剖面越來越接近于對(duì)數(shù)律,說明剖面逐漸向完全湍流發(fā)展。

4 結(jié) 論

本文以馬赫數(shù)為6的后掠鈍板邊界層為研究對(duì)象,采用直接數(shù)值模擬方法研究了橫流定常渦從首次失穩(wěn)、二次失穩(wěn)到轉(zhuǎn)捩發(fā)生的過程。得到以下結(jié)論:

1)橫流定常渦的非線性作用引起平均流修正,可使壁面摩擦系數(shù)曲線有一定程度的抬升,但是橫流渦飽和后壁面摩擦系數(shù)不再增大,轉(zhuǎn)捩不會(huì)發(fā)生。

2)在橫流渦非線性飽和(幅值約為0.2)的基礎(chǔ)上,發(fā)生二次失穩(wěn),產(chǎn)生高頻的不穩(wěn)定波。二次失穩(wěn)波由于非線性作用產(chǎn)生的低頻諧波及定常渦迅速增長(zhǎng),促使壁面摩擦系數(shù)急劇抬升,同時(shí)首次失穩(wěn)的飽和橫流渦結(jié)構(gòu)破碎,最終促使轉(zhuǎn)捩發(fā)生。

致謝：此項(xiàng)工作得到天津大學(xué)周恒院士的關(guān)心和指導(dǎo),天津大學(xué)趙磊博士提供了DNS程序及全局穩(wěn)定性分析程序,并在國(guó)家超算天津中心“天河一號(hào)”完成計(jì)算,特此致謝。