液氫空化流動特性研究

王國玉，陳泰然，黃 彪，潘俊屹

(北京理工大學(xué)機械與車輛學(xué)院，北京 100081)

液氫空化流動特性研究

王國玉，陳泰然，黃 彪，潘俊屹

(北京理工大學(xué)機械與車輛學(xué)院，北京 100081)

對液氫空化流動的實驗數(shù)據(jù)進行了統(tǒng)計分析，數(shù)值計算了大范圍溫度變化下液氫空化流動的變化特征。結(jié)果表明采用的數(shù)值計算方法計算結(jié)果與實驗數(shù)據(jù)吻合良好。隨著溫度變化，液氫空化流動存在兩種空化行為。當(dāng)溫度低于轉(zhuǎn)捩溫度時，空化行為的變化受液/汽密度比控制，空穴面積隨溫度的增加而增加；當(dāng)溫度高于轉(zhuǎn)捩溫度時，空化行為的變化熱力學(xué)效應(yīng)控制，空穴面積隨溫度的增加而減小?？栈瘮?shù)和來流速度越大，轉(zhuǎn)捩溫度越高。

液氫；空化流動；轉(zhuǎn)捩溫度

0 引言

航空航天事業(yè)的高速發(fā)展對運載工具的推力、功率密度和可靠性提出了越來越高的要求，液體火箭發(fā)動機正是當(dāng)前大運載的主力軍。由于發(fā)動機元件的設(shè)計受到尺寸和重量的限制，提高推力就需要盡可能提高渦輪泵的功率密度，這對渦輪泵的空化性能提出了嚴(yán)苛的要求，空化條件也一直是渦輪泵設(shè)計的主要限制條件之一[1]。常溫水體中的空化相變過程通常被假設(shè)為等溫過程[2]。然而，由于熱流體較小的液/汽相密度比和相變過程中顯著的溫度變化，等溫假設(shè)對于熱流體空化而言則不再成立[3]。液體火箭發(fā)動機所用的液氫和液氧等燃料就是典型的熱流體。研究表明[4]，誘導(dǎo)輪內(nèi)部液氫的非定?？栈鲃语@著影響誘導(dǎo)輪的性能，同時液氫空化的熱力學(xué)效應(yīng)影響顯著，從而改變了誘導(dǎo)輪內(nèi)部的流動特征。更加精確地描述強熱力學(xué)效應(yīng)作用下液體火箭發(fā)動機燃料的空化流動特性對液體火箭發(fā)動機的性能提升具有重要意義。

液氫和液氧等低溫流體的空化流動實驗研究[5]，對設(shè)備要求較高，花費較大，且危險系數(shù)高。目前只有NASA的研究人員Hord[6]在1973年做了最為完整的針對液氫和液氮的空化流動實驗，該實驗對認識低溫流體的空化流動特征提供了重要參考，也一直是發(fā)展數(shù)值計算方法的驗證指標(biāo)。近年來，隨著計算機水平的提高，空化流動的數(shù)值計算方法發(fā)展迅速?？栈Ｐ驮诳栈鲃拥臄?shù)值計算中扮演著重要的角色，計算中通常假設(shè)汽液混合相為均質(zhì)，然后利用正壓狀態(tài)方程模型或相間質(zhì)量傳輸方程模型來模擬兩相區(qū)域內(nèi)的變化。在正壓狀態(tài)方程模型中，混合物的密度由狀態(tài)方程確定，即混合區(qū)密度是壓力與溫度的函數(shù)。Rapposelli等[7]基于正壓狀態(tài)方程模型計算了液氮和液氧繞對稱回轉(zhuǎn)體和二維翼型的空化流場，利用正壓方程來計算空化流中的密度場，無法捕捉到斜壓矩的產(chǎn)生，尤其是沒有考慮相間質(zhì)量傳輸，無法解釋空化過程的對流和輸運，因而無法反映相變造成的熱力學(xué)效應(yīng)的影響機制。相間質(zhì)量傳輸方程模型根據(jù)蒸發(fā)和凝結(jié)源項來調(diào)節(jié)空化區(qū)域的質(zhì)量傳輸，可以真實反映由于相變造成的當(dāng)?shù)販囟鹊淖兓?。Hosangadi等[8-9]基于相間質(zhì)量傳輸模型，在全流場內(nèi)求解了能量方程，并耦合了流體介質(zhì)物質(zhì)屬性的變化，預(yù)測的繞翼型的壓強和溫度與實驗結(jié)果吻合良好。Utturkar等[1]采用基于界面捕捉空化模型來捕捉空化區(qū)域的變化，結(jié)果證明，翼型頭部溫度與實驗吻合較差。

空化流動是一種強烈的湍流流動，因此湍流模型的選擇對于模擬空化流動尤為重要?，F(xiàn)有的基于雷諾時均化(RANS)的湍流模型存在對空化區(qū)域湍流黏性過預(yù)測的問題[10-12]，無法準(zhǔn)確預(yù)測空化流動的非定常發(fā)展過程。近來發(fā)展的大渦模擬(LES)則可以很好地預(yù)測空化流動的非定常行為，但對計算資源要求較高[13-15]。因此，在現(xiàn)有湍流模型的基礎(chǔ)上發(fā)展了一系列混合模型，如基于密度修正的湍流模型(DCM)[11-12]、局部時均化湍流模型(PANS)[16]和濾波器湍流模型(FBM)[17]。

雖然國內(nèi)外研究人員已經(jīng)對液氫等低溫流體的空化流動特征開展了大量的實驗研究工作，但由于實驗條件限制，對于大范圍溫度變化條件下(尤其在三相點附近溫度)液氫的空化流動特性的變化還不十分清楚。

本文首先對Hord[6]完成的液氫繞軸對稱回轉(zhuǎn)體的空化流動實驗數(shù)據(jù)進行了統(tǒng)計處理，分析了在不同溫度、速度和空化數(shù)下空穴變化的特征。然后應(yīng)用基于傳熱修正的相間質(zhì)量傳輸模型和大渦模擬湍流模型，在求解連續(xù)方程、動量方程、質(zhì)量輸運方程的同時耦合能量方程，計算了液氫繞軸對稱回轉(zhuǎn)體的空化流動。通過與實驗數(shù)據(jù)對比驗證了數(shù)值計算方法的準(zhǔn)確性，然后計算了大范圍溫度變化條件下液氫空化流動的變化特征，為解決工程中誘導(dǎo)輪內(nèi)部低溫介質(zhì)的空化問題提供了重要參考。

1 控制方程及數(shù)值計算方法

1.1 控制方程

均質(zhì)多相流假設(shè)下空化流動的基本控制方程包括連續(xù)性方程、動量方程、能量方程和質(zhì)量方程，在笛卡兒坐標(biāo)下可以表示為：

(1)

(2)

(3)

(4)

ρm=ρl(1-αv)+ρvαv

(5)

μm=μl(1-αv)+μvαv

(6)

其中，ρm為混合密度,ρl和ρv分別為液相和汽相的密度,U為速度,p為壓力,μm為混合黏度,μl和μv分別為液相和汽相的動力黏度,μtur為湍流黏性系數(shù),f為質(zhì)量分?jǐn)?shù),Lev為蒸發(fā)潛熱,h為比焓,Prtur和Prlam分別為湍流和層流普朗特數(shù),αl和αv分別為液相和汽相的體積分?jǐn)?shù).

方程(4)中的m+和m-分別表示凝結(jié)和蒸發(fā)源項. 基于傳熱修正的相間質(zhì)量傳輸模型可以表示為[3]：

(7)

(8)

(9)

其中，M為摩爾質(zhì)量,Rv為氣體體積常數(shù),T為當(dāng)?shù)販囟?，pv為飽和蒸汽壓.C1= 0.13，C2= 0.01 為經(jīng)驗系數(shù)，該經(jīng)驗系數(shù)已在繞水翼和收縮擴張流道內(nèi)的空化流動計算中得到驗證[3,20-21]，因此本文直接采用。氣泡半徑rb根據(jù)單泡壁面?zhèn)鳠徇^程推導(dǎo)而來[3]，其中,Lev為蒸發(fā)潛熱,t∞為參考時間(根據(jù)特征長度和流動速度計算而來), ΔT為當(dāng)?shù)販囟茸兓?λl為液相熱擴散率,Cl為液相比熱容。

字際關(guān)系指字與字因來源或使用場景而產(chǎn)生的關(guān)系，包括古今字、異體字、繁簡字關(guān)系等。古代漢語與現(xiàn)代漢語中的字際關(guān)系不盡相同，同一個字記錄的詞語也有所不同。古代漢語中兩個音義不同的字，在現(xiàn)代漢語中可能“歸并”為同一個字。如“后”與“後”，在古代漢語中迥然有別，但在現(xiàn)代漢語中“後”被廢除，“后”一身而兼二義。整理古籍時，我們不能拿現(xiàn)代漢字的文字體系來對應(yīng)古代漢字的文字體系。特別是繁體字整理本，應(yīng)盡量保留古籍的用字原貌，不能對古籍用字妄加臆改，或者對字際關(guān)系視而不見。

1.2 大渦模擬湍流模型

本文采用大渦模擬對湍流進行處理，控制方程如下：

(10)

式中亞格子應(yīng)力τij定義如下：

τij=UiUj-UiUj

(11)

采用Smagorinsky 模型[18]對亞格子應(yīng)力進行處理

(12)

其中，Sij為大尺度應(yīng)變量張量，亞格子黏性νSGS通過以下方法處理

(13)

本文中計算域與邊界條件設(shè)置與Hord[6]實驗保持一致，考慮到計算的經(jīng)濟性，將軸對稱回轉(zhuǎn)體模型簡化為二維，計算域及邊界條件如圖1(a)所示。本文所采用的二維計算域存在一定誤差，但其可行性已在相關(guān)文獻中得到驗證，并不影響本文對液氫空化流動特性的分析[22-23]。為了更準(zhǔn)確地反映空化流動細節(jié)特征并滿足大渦模擬計算要求，如圖1(b)在回轉(zhuǎn)體近壁面周圍進行了網(wǎng)格加密，并確保近壁面第一層網(wǎng)格的y+值在1左右，滿足壁面函數(shù)的要求。邊界條件采用速度入口和壓力出口并在邊界處設(shè)置了溫度和相分布，流動區(qū)域上下邊界為自由滑移壁面條件，回轉(zhuǎn)體表面采用絕熱、無滑移壁面。

為了驗證計算網(wǎng)格，圖2給出了349B工況[6]下實驗數(shù)據(jù)與不同網(wǎng)格計算得到的溫度和壓力對比圖，邊界條件與Hord實驗[3]保持一致。從圖中可以看出，本文采用的數(shù)值計算方法能基本描述液氫繞對稱回轉(zhuǎn)體流動中溫度和壓力分布，與實驗數(shù)據(jù)基本一致。在空穴產(chǎn)生的區(qū)域內(nèi)，溫度和壓力降低，在空穴潰滅的區(qū)域，溫度略有升高，這與熱力學(xué)效應(yīng)對空化的影響效果是一致的。圖中顯示兩種網(wǎng)格數(shù)值計算結(jié)果與實驗均比較接近，因此在后續(xù)的計算中采用網(wǎng)格數(shù)較小的網(wǎng)格進行計算，網(wǎng)格總數(shù)為1.1×105。

2 結(jié)果與討論

2.1 液氫非定?？栈鲃犹卣?/p>

為直觀顯示液氫空化流動的非定常演化過程，圖3給出了考慮熱力學(xué)效應(yīng)和等溫條件下空穴發(fā)展過程，兩種條件下均以0.2tcycle(tcycle為空化發(fā)展周期)為時間間隔來展示動態(tài)演化過程。從圖中可以明顯看到，兩種條件下的空穴均為附著型空穴，均在回轉(zhuǎn)體前緣生長并逐漸向后發(fā)展，最后發(fā)生脫落。然而，在不同條件下，液氫非定?？栈鲃拥奶卣饕泊嬖诤艽蟛煌Ｔ诳紤]熱力學(xué)效應(yīng)條件下，空穴內(nèi)部大部分區(qū)域的氣體含量少于55%，最大能達到80%；而在等溫條件下，空穴內(nèi)大部分區(qū)域的氣體均達到了90%，最大為95%，且等溫條件下的空穴面積明顯更大。結(jié)合圖2(a)可以看出，在考慮熱力學(xué)效應(yīng)條件下，液氫蒸發(fā)吸熱使得當(dāng)?shù)販囟冉档?，從而改變了?dāng)?shù)匾簹涞奈镔|(zhì)屬性，當(dāng)?shù)仫柡驼羝麎航档鸵馕吨?dāng)?shù)乜栈瘮?shù)的升高。說明熱力學(xué)效應(yīng)抑制了液氫空化的發(fā)展，這是液氫空化較常溫水體空化最顯著的區(qū)別。進一步分析可以看出，在考慮熱力學(xué)效應(yīng)條件下，空穴內(nèi)外的溫度梯度和壓力梯度均明顯減小，空穴變得模糊，發(fā)展的尺度和程度均比等溫條件有所減小。顯然，耦合了能量方程的算法更能真實反映液氫空化過程，與實驗結(jié)果更為吻合。

2.2 溫度變化對液氫空化流動的影響

上文分析了熱力學(xué)效應(yīng)對液氫空化流動的顯著影響，接下來通過對Hord[6]完成的實驗數(shù)據(jù)進行統(tǒng)計處理，分析在不同溫度、速度和空化數(shù)條件下空穴變化的分布。圖4給出了基于實驗數(shù)據(jù)擬合得來的不同來流條件下的空穴長度分布，圖中黑點為實驗數(shù)據(jù)點。每幅圖中的數(shù)據(jù)均保證了相近的來流速度，從而得到在相同來流速度條件下不同空化數(shù)和溫度對空穴長度的影響。圖4(a)和圖4(b)中空穴面積的變化呈現(xiàn)相同的規(guī)律，在同一來流溫度下，空化數(shù)越大則空穴長度越??；在空化數(shù)相等的情況下，溫度越高則空穴長度越小。這與上文討論的熱力學(xué)效應(yīng)帶來的影響是相符的。隨著來流速度進一步提高，如圖4(c)和圖4(d)所示，空穴長度的整體變化基本保持不變。然而，在來流溫度為22.5K附近，出現(xiàn)了隨著溫度增加空穴長度不變甚至增加的區(qū)域，這就出現(xiàn)了所謂的“反熱力學(xué)效應(yīng)”區(qū)域。液氫空化流動中的“反熱力學(xué)效應(yīng)”現(xiàn)象與高溫水體中觀測得到的結(jié)果一致[3,19]，然而此前并未報道過。

為了進一步分析溫度變化對液氫空化流動的影響，本文通過數(shù)值計算得到了不同流動條件下的空化流動變化。圖5給出了空化發(fā)展10個準(zhǔn)周期內(nèi)無量綱空穴平均面積隨溫度的變化。從圖5(a)中來流速度U∞=63.9m/s的工況可以看出，在該工況下空穴面積隨著來流溫度增加單調(diào)減?。欢鴮τ趤砹魉俣萓∞=90m/s的工況，空穴面積隨著來流溫度增加出現(xiàn)了先增后減的變化過程。結(jié)合高溫水體的分析結(jié)果[3]和液氫的物質(zhì)屬性變化，認為當(dāng)液氫的溫度接近其三相點時，空化的熱力學(xué)效應(yīng)并不十分顯著，此時液/汽相的密度比主導(dǎo)了空化的行為。隨著溫度的增加，液/汽相密度比逐漸減小，空穴內(nèi)部的小氣泡直徑變小，空穴形態(tài)變得模糊，空穴內(nèi)部的小碎泡以更快的速度聚合并向下游移動，從整體上造成了空穴面積的增加。隨著溫度的進一步增加，液氫空化的熱力學(xué)效應(yīng)顯著增加，液相蒸發(fā)吸熱使得當(dāng)?shù)販囟冉档?，空穴區(qū)域內(nèi)溫度的降低使得液相的飽和蒸汽壓降低，飽和蒸汽壓降低意味空化更難發(fā)生，這就是空化熱力學(xué)效應(yīng)抑制空化發(fā)展的基本過程。因此，對于σ∞=0.38，U∞=90m/s的工況，空化行為轉(zhuǎn)變的轉(zhuǎn)捩溫度在16K左右。進一步分析空穴面積在大空化數(shù)下的變化(如圖5(b)和圖5(c)所示)，可以發(fā)現(xiàn)空化數(shù)越大，則轉(zhuǎn)捩溫度越高，說明空化數(shù)越大空化熱力學(xué)效應(yīng)越弱。對比不同速度下空穴的變化可以發(fā)現(xiàn)，來流速度增加使得空化流動慣性力增加，從而削弱了熱力學(xué)效應(yīng)對空穴的影響。因此，在相同空化數(shù)下，來流速度越大，轉(zhuǎn)捩溫度越高。

本文的研究工作為誘導(dǎo)輪的設(shè)計研制工作提供了重要參考。工程中應(yīng)該根據(jù)誘導(dǎo)輪的性能需求，合理設(shè)置液氫溫度、入口壓力、流量和誘導(dǎo)輪轉(zhuǎn)速的分配，避免誘導(dǎo)輪的運行工況在轉(zhuǎn)捩溫度附近運行。

3 結(jié)論

本文首先對液氫繞軸對稱回轉(zhuǎn)體的空化流動實驗數(shù)據(jù)進行了統(tǒng)計處理，分析了在不同溫度、速度和空化數(shù)下空穴變化的特征。然后應(yīng)用基于傳熱修正的相間質(zhì)量傳輸模型和大渦模擬湍流模型，在求解連續(xù)方程、動量方程、質(zhì)量輸運方程的同時耦合能量方程，計算了大范圍溫度變化條件下液氫空化流動的變化特征，本文主要結(jié)論有以下幾點：

1)本文采用的數(shù)值計算方法計算結(jié)果與實驗值吻合良好，能夠較準(zhǔn)確地預(yù)測液氫繞回轉(zhuǎn)體空化流動的非定常演變過程，準(zhǔn)確反映熱力學(xué)效應(yīng)對液氫空化流動的影響。

2)隨著溫度的變化，液氫空化流動存在兩種不同的空化行為，空化行為轉(zhuǎn)變所對應(yīng)的溫度為該工況下的轉(zhuǎn)捩溫度。當(dāng)溫度低于轉(zhuǎn)捩溫度時，空化行為的變化受液/汽相的密度比控制，空穴面積隨著溫度的增加而增加；當(dāng)溫度高于轉(zhuǎn)捩溫度時，空化行為的變化熱力學(xué)效應(yīng)控制，空穴面積隨著溫度的增加而減小。

3)在其他流動條件不變的前提下，空化數(shù)和來流速度越大，則轉(zhuǎn)捩溫度越高。

4)工程中應(yīng)根據(jù)誘導(dǎo)輪的性能需求，合理設(shè)置液氫溫度、入口壓力、流量和誘導(dǎo)輪轉(zhuǎn)速的分配，避免誘導(dǎo)輪在轉(zhuǎn)捩溫度附近運行。

[1] Utturkar Y, Wu J, Wang G, et al. Recent progress in modeling of cryogenic cavitation for liquid rocket propulsion[J].Progress in Astronautical Sciences, 2005, 41(7): 558-608.

[2] 張博,王國玉,黃彪,等.云狀空化非定常脫落機理的數(shù)值與實驗研究[J].力學(xué)學(xué)報,2009, 41(5):652-658. Zhang Bo, Wang Guo-yu, Huang Biao, et al. Numerical and experimental studies on unsteady shedding mechanisms of cloud cavitation[J].Chinese Journal of Theoretical and Applied Mechanics,2009,41(5):652-658(in Chinese).

[3] Chen T, Huang B, Wang G.. Numerical study of cavitating flows in a wide range of water temperatures with special emphasis on two typical cavitation dynamics[J].International Journal of Heat and Mass Transfer, 2016, 101: 886-900.

[4] Hosangadi A, Ahuja V, Ungewitter R J, et al. Analysis of thermal effects in cavitating liquid hydrogen inducers[J].Journal of Propulsion and Power, 2007, 23: 1225-1234.

[5] Franc J, Rebattet C, Coulon A. An experimental investigation of thermal effects in a cavitating inducer[J].Journal of Fluids Engineering,2004, 126: 716-723.

[6] Hord J. Cavitation in liquid cryogenics[R].III-Ogive, NASA-CR-2242,1973.

[7] Rapposelli E, d’Agostino L. A barotropic cavitation model with thermodynamic effects[C].Fifth International Symposium on Cavitation, Osaka, Japan, 2003: 1-9.

[8] Hosangadi A, Ahuja V. Numerical study of cavitation in cryogenic fluids[J].Journal of Fluid Engineering,2005,127:267-281.

[9] Hosangadi A, Ahuja V, Ungewitter R J, et al. Numerical study of a flat plate inducer: comparison of performance in liquid hydrogen and water[C].Proceedings of the 42ndAIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Sacramento, CA, United States,2006:7210-7222.

[10] Reboud J L, Stutz B, Coutier O. Two phase flow structure of cavitation:experiment and modeling of unsteady effects[C]. 3rdInternational Symposium on Cavitation (CAV1998), Grenoble, France, 1998.

[11] Coutier-Delgosha O, Fortes-Patella R, Reboud J L. Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation[J].Journal of Fluids Engineering,2003, 125 (1):38-45.

[12] Huang B, Dudoin A, Young Y L. Physical and numerical investigation of cavitating flows around a pitching hydrofoil[J].Physics of Fluids, 2013,25 (10): 102109.

[13] Wang G, Ostoja-Starzewski M. Large eddy simulation of a sheet/cloud cavitation on a NACA0015 hydrofoil[J].Applied Mathematical Modelling,2007, 31 (3):417-447.

[14] Ji B, Luo X W, Arndt R E A, et al. Large eddy simulation and theoretical investigations of the transient cavitating vortical flow structure around a NACA66 hydrofoil[J].International Journal of Multiphase Flow,2015,68: 121-134.

[15] Huang B, Zhao Y, Wang G. Large eddy simulation of turbulent vortex-cavitation interactions in transient sheet/cloud cavitating flows[J]. Computers & Fluids, 2014, 92 (3) :113-124.

[16] Ji B, Luo X, Wu Y, et al. Numerical analysis of unsteady cavitating turbulent flow and shedding horse-shoe vortex structure around a twisted hydrofoil[J]. International Journal of Multiphase Flow, 2013, 51(5):33-43.

[17] Wu J Y, Wang G Y, Shyy W. Time-dependent turbulent cavitating flow computations with interfacial transport and filter-based models[J].International Journal for Numerical Methods in Fluids,2005, 49 (7): 739-761.

[18] Smagorinsky J. General circulation experiments with the primitive equations[J].Monthly weather review, 1963, 91(3): 99-164.

[20] Liu S, Li S, Zhang L, et al. A mixture model with modified mass transfer expression for cavitating turbulent flow simulation[J].Engineering Computations,2008,25 (4): 290-304.

[21] Liu D M, Liu S H, Wu Y L, et al. A thermodynamic cavitation model applicable to high temperature flow[J].Thermal Science, 2011, 15(suppl. 1):95-101.

[22] Zhu J, Zhao D, Lu X,et al. Interactions of vortices, thermal effects and cavitation in liquid hydrogen cavitating flows[J].International Journal of Hydrogen Energy, 2016, 41(1):614-631

[23] Long X, Liu Q, Ji B. Numerical investigation of two typical cavitation shedding dynamics flow in liquid hydrogen with thermodynamic effects[J].International Journal of Heat and Mass Transfer,2017, 109: 879-893.

Characteristics of Hydrogen Cavitating Flows in a WideRange of Free-stream Temperature

WANG Guo-yu, CHEN Tai-ran, HUANG Biao, PAN Jun-yi

(School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China)

In this paper, the experimental data of liquid hydrogen cavitating flow are statistically analyzed, and the characteristics of liquid hydrogen cavitating flow under a wide range of temperatures are numerically investigated. The results show that the calculated results are in good agreement with the experimental data. As the temperature changes, there are two kinds of cavitation behavior in liquid hydrogen. The change of the cavitation behavior is dominated by the liquid/vapor density ratio when the temperature is lower than the transition temperature, and the cavity area increases with the increasing temperature. When the temperature is higher than the transition temperature, the thermodynamic effects dominate cavitation behavior, and the cavity area decreases with the increasing temperature. Moreover, the transition temperature increases with the increasing cavitation number and free-stream velocity.

Hydrogen; Cavitating flows; Transition temperature

符號表

2017-03-16；

2017-04-21基金項目：國家自然科學(xué)基金(51479002)

王國玉(1961-)，男，教授，博導(dǎo)，主要從事空化與水動力學(xué)研究。E-mail:wangguoyu@bit.edu.cn

V419

A

2096-4080(2017)01-0027-07