非設(shè)計(jì)工況下混流式水輪機(jī)尾水管壓力脈動(dòng)研究

張 廣,許 彬,魏顯著

(哈爾濱大電機(jī)研究所 水力發(fā)電設(shè)備國(guó)家重點(diǎn)試驗(yàn)室，哈爾濱 150040)

非設(shè)計(jì)工況下混流式水輪機(jī)尾水管壓力脈動(dòng)研究

張 廣,許 彬,魏顯著

(哈爾濱大電機(jī)研究所 水力發(fā)電設(shè)備國(guó)家重點(diǎn)試驗(yàn)室，哈爾濱 150040)

本文對(duì)非設(shè)計(jì)工況下水輪機(jī)尾水管流動(dòng)進(jìn)行數(shù)值研究。通過(guò)求解雷諾時(shí)均N-S方程和DES(Detached Eddy Simulation)湍流方程，分析了尾水渦帶的形態(tài)及其引起的壓力脈動(dòng)特性。研究結(jié)果表明，尾水渦帶的長(zhǎng)度、偏心距及轉(zhuǎn)動(dòng)方向與水輪機(jī)運(yùn)行條件有關(guān)：錐管上游測(cè)點(diǎn)周期性壓力脈動(dòng)比較規(guī)律，壓力脈動(dòng)幅值基本相同且主頻突出；錐管下游同一平面上的不同測(cè)點(diǎn)壓力脈動(dòng)的幅值差異較大，而且壓力脈動(dòng)的頻率成分增多。

混流式水輪機(jī)；非設(shè)計(jì)工況；壓力脈動(dòng)；渦帶

隨著電力市場(chǎng)需求的變化水電廠調(diào)峰任務(wù)凸顯，水輪機(jī)經(jīng)常需要在偏離設(shè)計(jì)工況條件下運(yùn)行，這時(shí)尾水管內(nèi)便會(huì)產(chǎn)生螺旋形的旋轉(zhuǎn)渦帶，渦帶所引起的壓力脈動(dòng)是造成機(jī)組振動(dòng)的主要原因，嚴(yán)重威脅著水電機(jī)組的安全運(yùn)行[1-2]。多年來(lái)，國(guó)內(nèi)外學(xué)者就有關(guān)尾水管渦帶引起的壓力脈動(dòng)問(wèn)題進(jìn)行了大量的研究。研究人員通過(guò)現(xiàn)場(chǎng)試驗(yàn)和模型試驗(yàn)手段，分析了不同工況下尾水管壓力脈動(dòng)信號(hào)的特征，對(duì)尾水管中渦帶形態(tài)的發(fā)展變化進(jìn)行觀測(cè)，探討了減小尾水管壓力脈動(dòng)的方法和手段[3-5]。但是由于尾水渦帶流動(dòng)現(xiàn)象的復(fù)雜性，試驗(yàn)研究在很多方面受到客觀條件的限制，難以獲得流場(chǎng)內(nèi)部物理量的分布情況。在數(shù)值仿真研究方面，Ciocan等人[6]使用RANS(雷諾平均)模型模擬了尾水管內(nèi)部流動(dòng)并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較，結(jié)果表明渦帶頻率和壓力脈動(dòng)幅值均存在偏差，同時(shí)越靠近轉(zhuǎn)動(dòng)中心軸向平均速度數(shù)值預(yù)測(cè)偏差越大。Vu.T.C等人[7]使用了k-ε兩方程模型數(shù)值模擬了尾水管內(nèi)部的流動(dòng)狀態(tài)和損失情況，研究發(fā)現(xiàn)水輪機(jī)運(yùn)行范圍一旦偏離最優(yōu)點(diǎn)，k-ε湍流模型的預(yù)測(cè)精度開始惡化，特別是在部分負(fù)荷條件下。Dhiman.S等人[8]測(cè)試了幾種RANS湍流模型在圓錐形擴(kuò)散段內(nèi)部的旋渦流動(dòng)的預(yù)測(cè)能力，結(jié)果表明URANS(非穩(wěn)態(tài)雷諾平均)模型均不能夠準(zhǔn)確的預(yù)測(cè)這種類型的流動(dòng)，特別是在距離進(jìn)口較遠(yuǎn)的尾水管下游。目前，LES(大渦模擬)方法也開始用于模擬尾水管渦帶流動(dòng)，LES方法在壓力脈動(dòng)頻率及幅度方面預(yù)測(cè)精度較高，但基于LES方法的非定常數(shù)值計(jì)算需要較為精細(xì)的網(wǎng)格和更小的時(shí)間步長(zhǎng)，對(duì)計(jì)算機(jī)性能的要求比較高[9]。數(shù)值仿真研究與試驗(yàn)研究可以在一定程度上彌補(bǔ)試驗(yàn)的不足，但是在數(shù)值預(yù)測(cè)精度及計(jì)算成本上有待于進(jìn)一步提高與平衡。

基于以上研究，本文使用DES(Detached Eddy Simulation)湍流模型對(duì)國(guó)內(nèi)某大型混流式水輪機(jī)模型機(jī)尾水管內(nèi)的非定常流動(dòng)進(jìn)行仿真。預(yù)測(cè)了不同運(yùn)行工況下尾水管內(nèi)部渦帶的形態(tài)及壓力脈動(dòng)情況，著重分析了尾水不同位置測(cè)點(diǎn)壓力脈動(dòng)信號(hào)的特點(diǎn)及規(guī)律。DES模型結(jié)合了RANS和LES模型的優(yōu)點(diǎn)，很好的平衡了數(shù)值預(yù)測(cè)精度及計(jì)算成本兩方面的矛盾，可以較好的預(yù)測(cè)非最優(yōu)設(shè)計(jì)工況條件下的非定常流場(chǎng)。

1 數(shù)值模型

1.1 基本控制方程

基本控制方程包括連續(xù)性方程和動(dòng)量方程如下：

(1)

(2)

式中:ρ、μ分別為按體積分?jǐn)?shù)確定的混合物的密度和動(dòng)力黏度;u為速度;S為源項(xiàng)。

1.2 湍流模型

本文數(shù)值仿真采用基于SST湍流模型的DES模型。DES模型結(jié)合了RANS模型和LES方法，其中RANS方法用于模擬邊界層內(nèi)以及輕微分離區(qū)域，LES用來(lái)模擬充分發(fā)展湍流區(qū)域。DES模型相對(duì)于LES方法放寬了對(duì)邊界層內(nèi)網(wǎng)格密度的要求，從而在高雷諾數(shù)流動(dòng)條件下極大地減少了網(wǎng)格數(shù)目，很大程度上節(jié)約了計(jì)算資源。

DES湍流模型在求解時(shí)，通過(guò)網(wǎng)格尺寸判定LES區(qū)域與RANS區(qū)域。當(dāng)SST方法預(yù)測(cè)湍流長(zhǎng)度比局部網(wǎng)格尺寸大時(shí)，將實(shí)現(xiàn)由SST模型向LES模型的轉(zhuǎn)換。這種情況下，在湍動(dòng)能輸運(yùn)方程中計(jì)算耗散率ε時(shí)采用的長(zhǎng)度由局部網(wǎng)格尺寸Δ取代，表達(dá)式如下：

ε=β*kω=k1.5/Lt

(3)

當(dāng)Lt>CDESΔ時(shí)

(4)

式中:k為湍動(dòng)能，ε為湍流耗散率，ω為湍動(dòng)頻率。

關(guān)于模型的詳細(xì)論述詳見(jiàn)參考文獻(xiàn)[10]。

1.3 計(jì)算對(duì)象及邊界條件

本文針對(duì)國(guó)內(nèi)某大型混流式水輪機(jī)模型進(jìn)行研究，建立了模型水輪機(jī)全通道三維計(jì)算模型，水輪機(jī)模型的基本參數(shù)見(jiàn)表1。

表1 水輪機(jī)模型基本參數(shù)Table 1 Basic parameters of water turbine model

水輪機(jī)模型包括蝸殼、導(dǎo)葉區(qū)、轉(zhuǎn)輪和尾水管四個(gè)部分，水輪機(jī)模型如圖1所。數(shù)值建模中將以上每個(gè)部分獨(dú)立的進(jìn)行網(wǎng)格劃分，之后通過(guò)交界面進(jìn)行連接，本文計(jì)算域網(wǎng)格單元總數(shù)約為865萬(wàn)。蝸殼及固定導(dǎo)葉部分采用四面體網(wǎng)格劃分，其它過(guò)流部件均采用六面體結(jié)構(gòu)化網(wǎng)格劃分，轉(zhuǎn)輪葉片部分的網(wǎng)格如圖2所示。

圖1 水輪機(jī)仿真模型Fig.1 Water turbine simulation model

圖2 轉(zhuǎn)輪葉片網(wǎng)格劃分Fig.2 Runner blade mesh generation

邊界條件方面，蝸殼進(jìn)口采用質(zhì)量流量進(jìn)口條件，尾水管出口采用壓力出口條件，根據(jù)模型試驗(yàn)中尾水管壓力初始化尾水管出口的壓力分布。動(dòng)域與靜域間的交界面(Interface)使用非一致網(wǎng)格連接，采用Transient Rotor Stator界面?zhèn)鬟f模型進(jìn)行模擬，其它邊界如蝸殼、轉(zhuǎn)輪、尾水管表面均采用無(wú)滑移壁面條件。數(shù)值計(jì)算模型在尾水管直錐段布置了4個(gè)壓力脈動(dòng)測(cè)點(diǎn)分別位于尾水錐管P1和P2平面上，如圖3所示。P1和P2平面距導(dǎo)葉中心線的距離分別為0.74D1和1.30D1，這主要是根據(jù)IEC 60193標(biāo)準(zhǔn)確定。本文數(shù)值計(jì)算工況點(diǎn)選擇為額定水頭下0.7QBEP和1.25QBEP(QBEP為最優(yōu)設(shè)計(jì)工況點(diǎn)對(duì)應(yīng)的流量)，這兩個(gè)流量下尾水渦帶的演變和形態(tài)很有代表性。非穩(wěn)態(tài)數(shù)值計(jì)算以穩(wěn)態(tài)數(shù)值計(jì)算結(jié)果為初始條件。數(shù)值計(jì)算總時(shí)間t=0.976 s，為轉(zhuǎn)輪旋轉(zhuǎn)15轉(zhuǎn)所需的時(shí)間；計(jì)算時(shí)間步長(zhǎng)Δt=0.001 s，約為轉(zhuǎn)輪每旋轉(zhuǎn)6°所需要的時(shí)間。

圖3 尾水管壓力測(cè)點(diǎn)布置圖Fig.3 Layout of draft tube pressure measuring points

2 計(jì)算結(jié)果與分析

2.1 尾水渦帶形態(tài)對(duì)比及壓力場(chǎng)分析

在非定常CFD計(jì)算結(jié)果的后處理中，旋轉(zhuǎn)渦帶形態(tài)采用等壓力面的方式進(jìn)行捕捉和表示。圖4給出了水輪機(jī)在0.7QBEP和1.25QBEP兩種工況下尾水管內(nèi)部的旋轉(zhuǎn)渦帶形態(tài)?？梢钥吹轿菜軠u帶起始于轉(zhuǎn)輪泄水錐的低部，渦帶在自身旋轉(zhuǎn)的同時(shí)隨著主流向下游旋進(jìn)，呈現(xiàn)一定的周期性運(yùn)動(dòng)特征。

圖4 尾水管渦帶形態(tài)Fig.4 Draft tube vortex shapes

對(duì)比兩種條件下渦帶形態(tài)可以明顯發(fā)現(xiàn)，在0.7QBEP工況下尾水渦帶的直徑和渦帶偏心量較大，渦帶一直向下游延伸至尾水管的彎肘段部分。而在1.25QBEP工況下尾水渦帶僅存在于尾水錐管部分，沒(méi)有進(jìn)一步向下游發(fā)展，但該工況下旋轉(zhuǎn)渦帶的螺旋特征較為明顯。

圖5給出在兩種運(yùn)行工況下一個(gè)渦帶旋轉(zhuǎn)周期內(nèi)P1、P2中間橫截面的壓力分布情況?？梢钥吹轿菜疁u帶區(qū)域?yàn)槊黠@的橢圓形低壓區(qū)，在0.7QBEP條件下渦帶低壓區(qū)的范圍較大，渦帶的偏心距較大，渦帶低壓區(qū)已經(jīng)接觸到尾水管壁面，這種情況會(huì)導(dǎo)致尾水管壁面產(chǎn)生較為明顯的振動(dòng)和噪聲。在1.25QBEP條件下低壓渦帶影響的范圍較小，同時(shí)渦帶的偏心量較小，低壓區(qū)沒(méi)有明顯的影響到尾水管壁面，尾水管內(nèi)壓力脈動(dòng)程度相對(duì)要小，這在下一節(jié)會(huì)詳細(xì)論述。

圖5 尾水管渦帶壓力云圖Fig.5 Stress nephogram of draft tube vortex

從兩種條件下尾水渦帶隨時(shí)間的運(yùn)動(dòng)情況來(lái)看，在0.7QBEP條件下尾水渦帶旋轉(zhuǎn)方向與轉(zhuǎn)輪轉(zhuǎn)動(dòng)的方向相同，而在1.25QBEP工況下尾水渦帶旋轉(zhuǎn)方向與轉(zhuǎn)輪轉(zhuǎn)動(dòng)方向相反。其原因可通過(guò)圖6轉(zhuǎn)輪葉片出口速度三角形加以解釋，圖中u為轉(zhuǎn)輪的圓周速度、w為水流的相對(duì)速度、v為水流絕對(duì)速度。水輪機(jī)在最優(yōu)設(shè)計(jì)工況下運(yùn)行時(shí)，轉(zhuǎn)輪出口水流基本為法向出流，如圖中v2所示。當(dāng)水輪機(jī)偏離最優(yōu)工況在0.7QBEP和1.25QBEP條件下運(yùn)行時(shí)，轉(zhuǎn)輪出口水流的絕對(duì)速度會(huì)如圖中v1和v3所示，這兩種條件下水流絕對(duì)速度的圓周分量相反，繼而使得渦帶呈現(xiàn)相反的轉(zhuǎn)動(dòng)方向。

圖6 轉(zhuǎn)輪葉片出口速度三角形Fig.6 Runner blade outlet velocity triangle

2.2 尾水管渦帶壓力脈動(dòng)特點(diǎn)分析

下面進(jìn)一步對(duì)尾水管不同測(cè)點(diǎn)的壓力脈動(dòng)信號(hào)進(jìn)行分析，從而揭示尾水管渦帶的壓力脈動(dòng)特性。圖7、圖8給出0.7QBEP條件下不同測(cè)點(diǎn)壓力隨時(shí)間的變化情況及FFT頻譜分析結(jié)果。

從圖7中可以看到，尾水管測(cè)點(diǎn)P11、P12的壓力脈動(dòng)信號(hào)呈現(xiàn)典型的正弦周期性變化，兩個(gè)測(cè)點(diǎn)的壓力脈動(dòng)幅值相當(dāng)，其數(shù)值約為2.8 kPa。此外由于兩個(gè)測(cè)點(diǎn)的相對(duì)位置的關(guān)系，使得旋轉(zhuǎn)渦帶低壓區(qū)位于P11測(cè)點(diǎn)的同時(shí)，P12測(cè)點(diǎn)為高壓區(qū)，壓力信號(hào)存在半個(gè)周期的相位差。壓力信號(hào)的頻率特性如圖8所示，可以看到兩個(gè)測(cè)點(diǎn)的壓力脈動(dòng)主頻突出，其數(shù)值均為f=3.61 Hz，該頻率也正是尾水管內(nèi)渦帶的轉(zhuǎn)動(dòng)頻率，約為轉(zhuǎn)頻fn=15.37 Hz的0.235倍，在壓力信號(hào)的頻域圖中并沒(méi)有明顯的諧頻出現(xiàn)。在位于尾水錐管下游測(cè)點(diǎn)P21、P22壓力信號(hào)同樣表現(xiàn)出明顯的周期性特征，兩個(gè)測(cè)點(diǎn)壓力脈動(dòng)存在半個(gè)周期的相位差。所不同的是P21、P22兩個(gè)測(cè)點(diǎn)的壓力脈動(dòng)幅值相差較大，P22測(cè)點(diǎn)壓力脈動(dòng)幅值與上游測(cè)點(diǎn)P11、P12相當(dāng)，而P21的壓力脈動(dòng)幅值約為其它測(cè)點(diǎn)的2倍。這說(shuō)明尾水渦帶引起的壓力脈動(dòng)在尾水錐管下游不同位置的脈動(dòng)幅值存在很大差異。從壓力脈動(dòng)的頻率特性上可以看到，P21、P22測(cè)點(diǎn)壓力脈動(dòng)主頻f=3.61 Hz，同時(shí)渦帶頻率的倍頻成分增多，其幅值大小隨渦帶頻率倍數(shù)的升高而降低。

圖7 壓力脈動(dòng)時(shí)域圖Fig.7 Time-domain diagram of pressure pulsation

圖8 壓力脈動(dòng)頻域圖Fig.8 Frequency-domain diagram of pressure pulsation

圖9給出1.25QBEP條件下不同測(cè)點(diǎn)壓力隨時(shí)間的變化情況。該工況下渦帶的偏心量及長(zhǎng)度較小，渦帶引起的壓力脈動(dòng)對(duì)壁面測(cè)點(diǎn)的影響較弱，可以看到該工況下尾水管壓力脈動(dòng)幅值整體較小。上游測(cè)點(diǎn)P11、P12的壓力脈動(dòng)幅值相當(dāng)，約為1.5 kPa。下游P21、P22測(cè)點(diǎn)壓力脈動(dòng)呈現(xiàn)出類似于0.7QBEP工況下的幅值特點(diǎn)，P21測(cè)點(diǎn)的脈動(dòng)幅值約為P22測(cè)點(diǎn)的2倍。另外有所不同的是從各個(gè)測(cè)點(diǎn)壓力數(shù)值上來(lái)看，位于尾水管左側(cè)的P11、P21測(cè)點(diǎn)的壓力值整體上要高于同一平面上右側(cè)的P12、P22測(cè)點(diǎn)。從圖10 FFT頻譜分析來(lái)看該工況下的壓力脈動(dòng)特征頻率成分增多，上游測(cè)點(diǎn)P11、P12主頻f=4.52 Hz。下游測(cè)點(diǎn)位于旋轉(zhuǎn)渦帶尾部，由于旋轉(zhuǎn)渦帶的耗散形成許多小尺度的湍流脈動(dòng)，該位置流動(dòng)比較復(fù)雜且具有一定的隨機(jī)性，這一點(diǎn)從頻率分析上也可以看到，在P21、P22測(cè)點(diǎn)有三個(gè)特征頻率分別為f1=4.52 Hz，f2=9.02 Hz，f3=14.52 Hz，而且信號(hào)中還夾雜著其它頻率成分。

圖9 壓力脈動(dòng)時(shí)域圖Fig.9 Time-domain diagram of pressure pulsation

圖10 壓力脈動(dòng)頻域圖Fig.10 Frequency-domain diagram of pressure pulsation

3 結(jié) 論

本文對(duì)混流式水輪機(jī)非設(shè)計(jì)工況下尾水渦帶形態(tài)及壓力脈動(dòng)進(jìn)行了數(shù)值研究，得到的主要結(jié)論如下：

1)水輪機(jī)尾水渦帶起始于轉(zhuǎn)輪泄水錐下端面，渦帶在自身旋轉(zhuǎn)的同時(shí)隨著主流向下游旋進(jìn)，尾水渦帶的長(zhǎng)度、偏心距及轉(zhuǎn)動(dòng)方向與水輪機(jī)運(yùn)行條件有關(guān)。

2)相對(duì)來(lái)說(shuō)，尾水錐管上游測(cè)點(diǎn)的壓力脈動(dòng)比較規(guī)律，壓力脈動(dòng)幅值相當(dāng)；錐管下游同一平面上的不同測(cè)點(diǎn)壓力脈動(dòng)的幅值存在較大差異，而且壓力脈動(dòng)特征頻率成分較多。

3)水輪機(jī)在部分負(fù)荷(0.7QBEP)條件下，尾水管內(nèi)低頻壓力脈動(dòng)幅值較大、主頻突出，水電機(jī)組應(yīng)盡量避免在該條件下運(yùn)行，必要時(shí)可采用外界擾動(dòng)的方法破壞尾水渦帶，減弱其壓力脈動(dòng)幅值和規(guī)律性，保證電站安全穩(wěn)定運(yùn)行。

[1] 黃源芳,劉光寧,樊世英.原型水輪機(jī)運(yùn)行研究[M].北京: 中國(guó)電力出版社,2010.HUANG Yuanfang,LIU Guangning,FAN Shiying.Research on operation of prototype turbine[M].Beijing: China Electric Power Press,2010.

[2] 吳玉林,吳曉晶,劉樹紅.水輪機(jī)內(nèi)部渦流與尾水管壓力脈動(dòng)相關(guān)性分析[J].水力發(fā)電學(xué)報(bào),2007,26(5):122-127.WU Yulin,WU Xiaojin,LIU Shuhong.Correlation analysis of the internal vortex flow of hydro turbine with the pressure fluctuation inside draft tube[J].Journal of Hydroelectric Engineering,2007,26(5): 122-127.

[3] 張飛,高忠信,潘羅平,等.混流式水輪機(jī)部分負(fù)荷下尾水管壓力脈動(dòng)試驗(yàn)研究[J].水利學(xué)報(bào),2011,42(10):1234-1238.ZHANG Fei,GAO Zhongxin,PAN Luoping,et al.Study on pressure fluctuation in Francis turbine draft tubes during partial load[J].Journal of Hydraulic Engineering,2011,42(10): 1234-1238.

[4] 何成連,王正偉,邱華.水輪機(jī)尾水管內(nèi)部壓力脈動(dòng)試驗(yàn)研究[J].機(jī)械工程學(xué)報(bào),2002,38(11):62-65.

HE Chenglian,WANG Zhengwei,QIU Hua.Experimental study on pressure surge in hydraulic turbine draft tube[J].Chinese Journal of Mechanical Engineering,2002,38(11): 62-65.

[5] Alin Ilie Bosioc,Romeo Susan-Resiga et al.Unsteady Pressure Analysis of a Swirling Flow With Vortex Rope and Axial Water Injection in a Discharge Cone.[J].Journal of Fluids Engineering,2012,134: 1-11.

[6] Ciocan,G.D.,Iliescu,M.S.,Vu,T.C.,Nennemann,B.,and Avellan,F.,2007,“Experimental Study and Numerical Simulation of the FLINDT Draft Tube Rotating Vortex,” ASME J.Fluids Eng.,129(2),pp.146-158.

[7] Vu,T.C.,Devals,C.,Zhang,Y.,Nennemann,B.,Guibault,F.,2011,“Steady and Unsteady Flow Computation in an Elbow Draft Tube With Experimental Validation,” Int.J.Fluid Mach.Syst.,4(1),pp.85-96.

[8] Dhiman.S.,Foroutan,H.,and Yavuzkurt,S.,2011,“Simulation of Flow Through Conical Diffusers With and Without Inlet Swirl Using CFD,” ASMEJSME-KSME Joint Fluids Engineering Conference,Hamamatsu,Japan,July 24-29,ASME Paper No.AJK2011-03005.

[9] Dragica Jost,Andrej Lipej.Numerical Prediction of Non-Cavitating and CavitatingVortex Rope in a Francis Turbine Draft Tube.Journal of Mechanical Engineering 2011,(57)6,445-456.

[10] A.Travin,M.Shur et al..Detached-eddy Simulations Past a Circular Cylinder .Internal Journal of Flow Turbulence and Combustion,2000,63:293-313.

(編輯 陳銀娥)

Study on pressure fluctuation in Francis turbine draft tubes under off-design condition

ZHANG Guang,XU Bin,WEI Xianzhu

(Harbin Institute of Large Electrical Machinery,State Key Laboratory of Hydropower Equipment,Harbin 150040,China)

In this paper,numerical study on the flow of water turbine draft tube in off-design conditions is carried out.The shape of the draft tube vortex and pressure pulsation characteristics induced by the vortex are analyzed by solving the Renault N-S equation and the DES (Detached Eddy Simulation) turbulence equation.The results show that the length,eccentricity and rotation direction of the draft tube vortex are related to the operating conditions of the turbine; the periodic pressure pulsation of the upper taper tube is relatively regulatory,amplitude of the pressure pulsation is basically the same and the main frequency is prominent.While,the amplitude of the pressure fluctuation at the same level of the lower taper tube is rather different,and the frequency component of the pressure pulsation increases.

Francis turbine; off-design working condition; pressure pulsation; vortex

2016-11-21;

2017-05-16。

張廣(1983—)，男，工程師，博士，主要從事水力機(jī)械優(yōu)化設(shè)計(jì)及多相流動(dòng)研究工作。

TK733+1

A

2095-6843(2017)04-0352-05