等離子體激勵(lì)位置對(duì)抑制壓氣機(jī)角區(qū)分離效果的影響

梁斐杰，陸利蓬，柳陽(yáng)威，孫槿靜

（北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京100191）

等離子體激勵(lì)位置對(duì)抑制壓氣機(jī)角區(qū)分離效果的影響

梁斐杰，陸利蓬，柳陽(yáng)威，孫槿靜

（北京航空航天大學(xué)能源與動(dòng)力工程學(xué)院，北京100191）

為了揭示等離子體氣動(dòng)激勵(lì)對(duì)角區(qū)分離的作用效果，應(yīng)用FLU EN T軟件數(shù)值模擬了等離子體激勵(lì)器對(duì)壓氣機(jī)葉柵角區(qū)分離的影響。采用等離子體激勵(lì)器的簡(jiǎn)化唯象模型，在壓氣機(jī)葉片吸力面和端壁不同位置沿流向施加激勵(lì)，對(duì)總壓損失系數(shù)、極限流線、不同截面流動(dòng)情況進(jìn)行了比較分析。結(jié)果表明：吸力面激勵(lì)對(duì)角區(qū)分離改善有限，角區(qū)未失速時(shí)，近分離點(diǎn)前是激勵(lì)最佳位置，角區(qū)失速后，激勵(lì)位置越靠前效果越好；端壁流向激勵(lì)能明顯減小角區(qū)分離損失，分離點(diǎn)至葉片前緣任何位置施加激勵(lì)效果一樣；組合激勵(lì)同時(shí)減小吸力面邊界層和端壁邊界層損失，使角區(qū)分離消失且不受攻角變化影響。

角區(qū)分離；流動(dòng)控制；等離子體激勵(lì)；唯象模型；數(shù)值模擬

0 引言

未來(lái)壓氣機(jī)發(fā)展的主要趨勢(shì)是更高的增壓比，更高的效率和更寬的穩(wěn)定工作范圍，但流動(dòng)分離卻導(dǎo)致了壓氣機(jī)性能的降低[1]。其中，在吸力面和端壁之間形成的角區(qū)分離流動(dòng)，是壓氣機(jī)內(nèi)部普遍存在的1種流動(dòng)分離現(xiàn)象，是壓氣機(jī)內(nèi)流動(dòng)損失和流動(dòng)堵塞的主要來(lái)源[2]，嚴(yán)重時(shí)引起失速和喘振[3]。因此，對(duì)壓氣機(jī)3維角區(qū)分離流動(dòng)機(jī)理、預(yù)測(cè)和控制的研究，始終是高性能先進(jìn)壓氣機(jī)設(shè)計(jì)中關(guān)注的重點(diǎn)問(wèn)題，也是內(nèi)流湍流研究的關(guān)鍵科學(xué)問(wèn)題[4-5]。

壓氣機(jī)角區(qū)分離是典型的3維分離，不同于2維分離，角區(qū)分離是在端壁及吸力面表面均有回流的區(qū)域，結(jié)構(gòu)復(fù)雜[6]。端壁出現(xiàn)明顯的大回流區(qū)時(shí)，又稱之為角區(qū)失速，會(huì)對(duì)性能造成更大的影響[7]，更需要進(jìn)行有效控制。近年來(lái)，流動(dòng)控制技術(shù)得到了快速發(fā)展，包括主動(dòng)控制技術(shù)和被動(dòng)控制技術(shù)2類(lèi)。其中，等離子體流動(dòng)控制技術(shù)是1種基于等離子體氣動(dòng)激勵(lì)的主動(dòng)控制方法，具有響應(yīng)迅速、作用頻帶寬、無(wú)需移動(dòng)機(jī)械部件、便于實(shí)時(shí)控制、功率消耗低等優(yōu)點(diǎn)。等離子體流動(dòng)控制技術(shù)在抑制翼型前緣流動(dòng)分離[8]、抑制渦輪[9-10]和壓氣機(jī)[11-13]流動(dòng)分離等方面已經(jīng)取得了一定的進(jìn)展。但相關(guān)控制機(jī)理和規(guī)律還有待進(jìn)一步研究，目前還沒(méi)有工程應(yīng)用中的壓氣機(jī)采用了等離子體激勵(lì)控制流動(dòng)分離。

本文采用FLUENT軟件數(shù)值模擬研究了不同攻角下等離子體激勵(lì)位置對(duì)抑制壓氣機(jī)葉柵角區(qū)分離效果的影響。數(shù)值計(jì)算中，等離子體激勵(lì)的模型采用Shyy等提出的唯象模型，通過(guò)FLUENT軟件自帶的UDF二次開(kāi)發(fā)模塊實(shí)現(xiàn)添加。研究算例采用劍橋大學(xué)的高負(fù)荷PVD葉柵，通過(guò)在葉片吸力面和端壁不同位置施加激勵(lì)，系統(tǒng)研究了等離子體氣動(dòng)激勵(lì)位置抑制壓氣機(jī)葉柵角區(qū)分離流動(dòng)的作用效果。

1 等離子體氣動(dòng)激勵(lì)原理及模型

1.1 等離子體氣動(dòng)激勵(lì)原理

等離子激勵(lì)器結(jié)構(gòu)[14]如圖1所示，將2層電極用絕緣介質(zhì)薄膜隔開(kāi)，其中1層電極完全暴露在空氣中，另外1層電極則被絕緣介質(zhì)所覆蓋。在2電極之間施加較高電壓時(shí)，電極間的氣體會(huì)發(fā)生電離，形成了大氣壓下均勻輝光放電等離子體，產(chǎn)生的等離子體對(duì)電極附近的帶電粒子作用體積力，從而使這些帶電粒子和中性的空氣粒子相互碰撞，并將動(dòng)量傳遞給空氣粒子，使其發(fā)生運(yùn)動(dòng)，最終對(duì)流動(dòng)產(chǎn)生加速作用。

圖1 等離子激勵(lì)器結(jié)構(gòu)

1.2 等離子體氣動(dòng)激勵(lì)數(shù)學(xué)模型

采用由Shyy等提出的唯象模型，該唯象模型[15]基于試驗(yàn)和理論分析，假設(shè)等離子體的作用區(qū)域如圖2所示，其中寬度為b，高度為a的直角三角形AOB，

圖2 等離子激勵(lì)器作用區(qū)域

式中A，B，C3項(xiàng)與通常的控制方程一樣，D項(xiàng)通過(guò)式（5）求解，即在該區(qū)域內(nèi)電場(chǎng)力呈線性分布，在三角形AOB以外的區(qū)域，電場(chǎng)強(qiáng)度的大小不足以擊穿空氣，所以三角形AOB區(qū)域?yàn)榈入x子體激勵(lì)器的有效擊穿區(qū)域。在原點(diǎn)O具有最大電場(chǎng)強(qiáng)度，在邊界OA、OB及AB上截?cái)嚯妶?chǎng)強(qiáng)度為Eb。

整個(gè)三角形區(qū)域里面電場(chǎng)力平行于AB邊界并呈線性變化，因此可以給出電場(chǎng)分布為

式中：U0為激勵(lì)電壓；d為電極間的間距；根據(jù)電場(chǎng)線性分布可知

通過(guò)式（1）中的電場(chǎng)強(qiáng)度分布可得該區(qū)域內(nèi)的電場(chǎng)力分布為

式中：?為施加電壓的頻率；α為電荷碰撞效率因子；ρc為電荷密度；ec為電子電荷常數(shù)；△t為激勵(lì)電壓的半周期；δ為狄拉克函數(shù)；Eb為電場(chǎng)邊界截止電壓。

為了模擬等離子體激勵(lì)器誘導(dǎo)的體積力對(duì)于流動(dòng)的作用效果，在N-S方程中添加相應(yīng)的體積力項(xiàng)。假設(shè)流動(dòng)是定常不可壓的低雷諾數(shù)流動(dòng)?？刂品匠虨?/p>

2 數(shù)值方法

2.1 算例簡(jiǎn)介

研究所選用的壓氣機(jī)葉柵是英國(guó)劍橋大學(xué)Gbadebo等人針對(duì)角區(qū)分離流動(dòng)機(jī)理研究所設(shè)計(jì)的高負(fù)荷平面葉柵，運(yùn)用試驗(yàn)和數(shù)值模擬方法系統(tǒng)研究了其內(nèi)部的3維角區(qū)分離流動(dòng)。該葉柵是1種典型的PVD（Prescribed Velocity Distribtion）葉柵，葉型主要參數(shù)見(jiàn)表1。

表1 葉柵信息

2.2 計(jì)算網(wǎng)格

計(jì)算網(wǎng)格由NUMECA軟件的AutoGrid模塊生成，采用分區(qū)結(jié)構(gòu)化網(wǎng)格，為O-4H型網(wǎng)格，網(wǎng)格塊間所有連接網(wǎng)格都是完全匹配的，以避免計(jì)算中引入插值誤差，所有壁面第1條網(wǎng)格線y+<1.4。由于該流動(dòng)具有對(duì)稱性，為減小計(jì)算量，只計(jì)算半個(gè)葉高，在葉中采用對(duì)稱邊界條件即可。為了減小網(wǎng)格對(duì)計(jì)算的影響，計(jì)算過(guò)程中數(shù)值試驗(yàn)了一系列網(wǎng)格，通過(guò)調(diào)整網(wǎng)格數(shù)和網(wǎng)格分布情況，得到網(wǎng)格無(wú)關(guān)解。最終使用的網(wǎng)格如圖3所示，展向網(wǎng)格數(shù)為51，網(wǎng)格總數(shù)944800。

圖3 計(jì)算網(wǎng)格

2.3 計(jì)算設(shè)置

計(jì)算采用FLUENT軟件，湍流模型選用RSM模型，壓力速度耦合求解采用SIMPLE算法，對(duì)流項(xiàng)和擴(kuò)散項(xiàng)分別采用2階迎風(fēng)格式和中心差分格式進(jìn)行離散。進(jìn)口給定速度分布，出口給定靜壓，葉中采用對(duì)稱邊界條件，固壁采用絕熱無(wú)滑移邊界條件。

2.4 等離子激勵(lì)器說(shuō)明

研究中等離子激勵(lì)器的作用區(qū)域?yàn)閍=1.5 mm，b=3 mm。施加電壓的頻率?=3 kHz；電荷密度ρc=1× 1011cm-3；電子電荷常數(shù)ec=1.602×10-19；放電時(shí)間△t=67 μs；電場(chǎng)邊界截止電壓Eb=30 kV/cm。

對(duì)等離子體激勵(lì)的數(shù)值模擬，通過(guò)FLUENT軟件的UDF（User Defined Functions）功能進(jìn)行二次開(kāi)發(fā)實(shí)現(xiàn)。等離子激勵(lì)器在不同位置施加方案如圖4所示，分別在葉片吸力面沿弦長(zhǎng)0%，30%，60%，90%及端壁對(duì)應(yīng)位置模擬等離子體氣動(dòng)激勵(lì)，從左到右分別為位置1～4。

圖4 等離子激勵(lì)器作用位置

3 結(jié)果及討論

3.1 不加等離子體氣動(dòng)激勵(lì)的流動(dòng)分析

為保證計(jì)算結(jié)果的可靠性，首先，對(duì)不加等離子體氣動(dòng)激勵(lì)0°攻角下的PVD葉柵進(jìn)行了校驗(yàn)，其中試驗(yàn)結(jié)果由英國(guó)劍橋大學(xué)Gbadebo等人完成[4]。端壁和吸力面的極限流線對(duì)比情況如圖5所示，從定性上看計(jì)算得到的角區(qū)分離和試驗(yàn)吻合得很好。

在不加等離子體氣動(dòng)激勵(lì)的情況下，對(duì)攻角為1°，3°，5°，7°的工況進(jìn)行了計(jì)算，并進(jìn)行極限流線對(duì)比如圖6所示。隨著攻角的增大，吸力面分離逐漸增大，3°攻角時(shí)，端壁出現(xiàn)很大的分離流動(dòng)，形成角區(qū)失速，使其性能急劇惡化。

3.2 吸力面不同位置施加激勵(lì)的流動(dòng)分析

圖5 0°攻角計(jì)算與試驗(yàn)極限流線對(duì)比

圖6 不同攻角下的極限流線對(duì)比

對(duì)應(yīng)每1個(gè)攻角的原始流動(dòng)，在吸力面的4個(gè)位置施加等離子體氣動(dòng)激勵(lì)，選擇攻角為1°和5°的工況進(jìn)行分析分別如圖7、8所示。對(duì)比發(fā)現(xiàn)，在吸力面施加激勵(lì)可以顯著改善吸力面上的分離流動(dòng)，使回流區(qū)面積減小，但對(duì)端壁的分離流動(dòng)改善有限。

首先，攻角為1°時(shí)，未發(fā)生角區(qū)失速，在位置3施加激勵(lì)可使吸力面分離幾乎消失；其次是在位置1和位置2施加激勵(lì)，二者效果相當(dāng)；再次，在位置4施加激勵(lì)效果不明顯。此時(shí)位置3更靠近分離點(diǎn)。當(dāng)攻角為5°時(shí)，角區(qū)發(fā)生失速，此時(shí)激勵(lì)位置越靠前效果越好。

激勵(lì)位置相同，攻角不同時(shí)，在位置1和位置2施加激勵(lì)都能使吸力面的流動(dòng)分離得到有效控制，端壁變化則不明顯。隨著激勵(lì)位置后移，激勵(lì)效果逐漸減弱，極限流線幾乎不變。

圖7 1°攻角時(shí)吸力面不同位置施加激勵(lì)極限流線對(duì)比

圖8 5°攻角時(shí)吸力面不同位置施加激勵(lì)極限流線對(duì)比

對(duì)比5°攻角，吸力面不同位置施加激勵(lì)時(shí)，葉片尾緣后50%弦長(zhǎng)位置截面的總壓損失，如圖9所示。

P01，P0，ρ，V1分別為進(jìn)口總壓，當(dāng)?shù)乜倝海べ|(zhì)密度和進(jìn)口速度。可見(jiàn)在吸力面不同位置施加激勵(lì)都不能使角區(qū)分離消失。激勵(lì)位置越靠前，吸力面邊界層的損失越小，角區(qū)中心的總壓損失系數(shù)越小。在前緣施加激勵(lì)能更好的抑制吸力面邊界層損失對(duì)角區(qū)分離的貢獻(xiàn)。

對(duì)比了不同攻角，不同位置施加激勵(lì)的總壓損失系數(shù)如圖10所示。計(jì)算同式（11）。攻角為0°和1°時(shí)未發(fā)生角區(qū)失速，此時(shí)位置1、2、3都在分離點(diǎn)附近或上游，施加激勵(lì)總壓損失系數(shù)變化基本一致。攻角增大到3°，發(fā)生角區(qū)失速，總壓損失急劇上升，激勵(lì)位置越靠近前緣，總壓損失系數(shù)越小，位置1施加激勵(lì)對(duì)不同攻角工況都能有效控制。

圖9 5°攻角時(shí)吸力面不同位置施加激勵(lì)時(shí)葉片尾緣50%弦長(zhǎng)截面的總壓損失

3.3 端壁不同位置施加激勵(lì)的流動(dòng)分析

對(duì)應(yīng)每1個(gè)攻角的原始流動(dòng)，在端壁的4個(gè)位置施加等離子體氣動(dòng)激勵(lì)，并選擇攻角為1°和5°的工況進(jìn)行分析分別如圖11、12所示。通過(guò)對(duì)比發(fā)現(xiàn)，端壁施加激勵(lì)可以顯著改善端壁上的分離流動(dòng)，甚至使分離消失。端壁流動(dòng)不再是由壓力面流向吸力面，而是在施加激勵(lì)的位置之后，流動(dòng)與主流相似。

圖10 總壓損失系數(shù)對(duì)比

圖111 °攻角時(shí)端壁不同位置施加激勵(lì)極限流線對(duì)比

圖12 5°攻角時(shí)端壁不同位置施加激勵(lì)極限流線對(duì)比

激勵(lì)位置越靠前，端壁流動(dòng)改善越好，位置1和2效果相似，同時(shí)會(huì)在吸力面形成回流區(qū)并隨攻角增大而增大。在靠后的位置施加激勵(lì)，隨著攻角增大，等離子體氣動(dòng)激勵(lì)對(duì)端壁和吸力面的流動(dòng)改善效果都變得有限，甚至并不會(huì)影響吸力面的流動(dòng)情況。

對(duì)比5°攻角，端壁不同位置施加激勵(lì)時(shí)，葉片尾緣后50%弦長(zhǎng)位置截面的總壓損失如圖13所示。發(fā)現(xiàn)在端壁位置1和2施加激勵(lì)，角區(qū)失速消失，端壁邊界層和吸力面邊界層沒(méi)有出現(xiàn)摻混，而是較好的保持了各自的獨(dú)立性，形成了損失相對(duì)較小的“2維”分離。在位置3、4施加激勵(lì)，角區(qū)失速并不會(huì)消失，尤其是在位置4施加激勵(lì)，變化有限?？梢?jiàn)激勵(lì)位置越靠前，激勵(lì)效果越好。

圖13 5°攻角時(shí)端壁不同位置施加激勵(lì)時(shí)葉片尾緣50%弦長(zhǎng)截面的總壓損失

在不同攻角下，不同位置施加激勵(lì)的總壓損失系數(shù)對(duì)比如圖14所示。發(fā)現(xiàn)在端壁施加激勵(lì)時(shí)，位置1和位置2效果幾乎一樣，但是隨著攻角的增大，位置1有優(yōu)于位置2的趨勢(shì)，這是因?yàn)榉蛛x點(diǎn)會(huì)隨攻角增大而前移。激勵(lì)位置越靠后，激勵(lì)效果越弱。同時(shí)，隨著攻角的增大，位置1和位置2的激勵(lì)效果急劇減弱，激勵(lì)效果對(duì)位置的敏感性變差。參考圖14說(shuō)明在大攻角工況下，吸力面邊界層2維分離帶來(lái)的損失加大。

圖14 總壓損失系數(shù)對(duì)比

3.4 吸力面和端壁的組合激勵(lì)

基于以上分析，在端壁位置1和吸力面位置1同時(shí)施加激勵(lì)。并對(duì)7°攻角，葉片尾緣后50%弦長(zhǎng)位置的周向平均總壓損失系數(shù)進(jìn)行了對(duì)比如圖15所示。采用組合激勵(lì)的方式，不論在近壁區(qū)還是主流區(qū)都可以產(chǎn)生很好的流動(dòng)控制效果，其結(jié)合了吸力面激勵(lì)和端壁激勵(lì)的優(yōu)點(diǎn)，使得在整個(gè)葉高范圍內(nèi)都有很好的作用效果。

圖15 7°攻角工況下周向平均總壓損失系數(shù)

4 結(jié)論

本文采用FLUENT軟件數(shù)值模擬研究了等離子體激勵(lì)位置對(duì)抑制高負(fù)荷壓氣機(jī)PVD葉柵角區(qū)分離效果的影響，得到主要結(jié)論如下：

（1）吸力面激勵(lì)可以改善吸力面邊界層并減小其對(duì)角區(qū)分離的貢獻(xiàn)，但對(duì)角區(qū)分離抑制效果有限；端壁激勵(lì)可以改善端壁邊界層流動(dòng)，可以較好抑制角區(qū)分離。

（2）未發(fā)生角區(qū)失速，近分離點(diǎn)施加激勵(lì)效果最優(yōu)；發(fā)生角區(qū)失速后，吸力面越靠前施加激勵(lì)抑制效果最好，端壁分離點(diǎn)至葉片前緣任何位置施加激勵(lì)抑制效果都基本一樣。

（3）端壁施加激勵(lì)對(duì)角區(qū)分離的抑制效果大于吸力面施加激勵(lì)，端壁邊界層對(duì)角區(qū)分離貢獻(xiàn)更大，是研究控制角區(qū)分離的重點(diǎn)。

（4）采取組合激勵(lì)的方式，可以很好抑制角區(qū)分離/失速，且不受攻角變化影響，可以據(jù)此發(fā)展更經(jīng)濟(jì)的組合方式進(jìn)行控制。

[1]Li Yinghong,Wu Yun,Zhou Min,et al.Control of the corner separation in a compressor cascade by steady and unsteady plasma aerodynamic actuation[J].Expro Fluids,2010,48: 1015-1023.

[2]Denton J D.Loss mechanisms in turbomachines[J].Journal of Turbomachinery,1993,115：621-656.

[3]Thiam A H,Whittlesey R W.Corner separation and the onset of stall in an axial compressor[R].AIAA-2008-4299.

[4]Gbadebo S A.Three-dimensional separations in compressors [D].United Kingdom:University of Cambridge,2003.

[5]Liu Yangwei,Lu Lipeng,Fang Le,et al.Modification of spalart-allmaras modelwith consideration ofturbulence energy backscatter using velocity helicity[J].Physics Letters A, 2011,375（24）：2377-2381.

[6]Joslyn H D,Dring R P.Axial compressor stator aerodynamics [J].Journal of Engineering for Gas Turbines and Power,1985, 107：485-493.

[7]Lei V M,Spakovszky Z S,Greitzer E M.A criterion for axial compressorhub-cornerstall[J].ASMEJournalofTurbomachinery, 2008,130（1）：1-10.

[8]Post ML,Corke T C.Separation control in high angle of attack airfoil using plasma actuators[J].AIAA Journal,2004,42（11）: 2177-2184.

[9]Huang J H,Corke T C.Plasma actuators for separation control of low pressure turbine blades[J].AIAA Journal,2006,44（1）：51-57.

[10]Huang J H,Corke T C.Unsteady plasma actuators for separation control of low-pressure turbine blades[J].AIAA Journal,2006,44（7）：1477-1487.

[11]吳云,李應(yīng)紅.等離子體氣動(dòng)激勵(lì)抑制壓氣機(jī)葉柵角區(qū)流動(dòng)分離的仿真與實(shí)驗(yàn) [J].航空動(dòng)力學(xué)報(bào),2009,24（4）：830-835. WU Yun,LI Yinghong.Computational and experimental of suppressingplasma aerodynamic actuation cornerflow separation control in a compressor cascade[J].Journal of Aerospace Power,2009,24（4）：830-835.(in Chinese)

[12]劉華坪，陳煥龍，袁繼來(lái)，等.變沖角條件下等離子體對(duì)擴(kuò)壓葉柵性能的影響 [J].航空動(dòng)力學(xué)報(bào),2011,26（8）：1871-1878. LIU Huaping,CHEN Huanlong,YUE Jilai,et al.Influence on the performance of compressor cascade with alterative incidence angles by the plasma[J].Journal of Aerospace Power,2011,26（8）：1871-1878.(in Chinese)

[13]趙小虎,吳云,李應(yīng)紅,等.高負(fù)荷壓氣機(jī)葉柵分離結(jié)構(gòu)及其等離子體流動(dòng)控制 [J].航空學(xué)報(bào),2012,33（2）: 208-219. ZHAO Xiaohu,WU Yun,LI Yinghong,et al.Separation structure and plasma flow controlon highly loaded compressure cascade[J].Journal of Aerospace,2012,33（2）: 208-219.(in Chinese)

[14]李綱.等離子體流動(dòng)控制機(jī)理極其應(yīng)用研究[D].北京：中國(guó)科學(xué)院工程熱物理研究所，2008. LI Gang.Investigation of plasma flow control mechanism and its application[D].Beijing:Chinese Academy of Sciences, 2008.(in Chinese)

[15]Shyy W,Jayaraman B.Modeling of glow dischargeinduced fluid dynamics[J].Journal of Applied Physics,2002,92（11）：6434-6443.

Impact of Plasma Actuating Position on Control of Corner Separation of a Compressor Cascade

LIANG Fei-jie，LU Li-peng，LIU Yang-wei，SUN Jin-jing
（School of Jet Propulsion,Beihang University,Beijing 100191,China）

To explore the intrinsic mechanism of plasma aerodynamic actuation on corner separation,FLUENT was used to simulate the effect of a plasma actuator on the corner separation of a compressor cascade.The phenomenological model was applied to impose actuation on different positions along with the flow direction both at the suction surface and hub.Total pressure loss coefficient,streamline,flow in different cross sections were comparatively analyzed.The results indicate that suction surface actuation reduces corner separation loss limited,the separation point is at the best actuating point without corner stall while the nearer actuation position is,the more effective is with corner stall.Hub actuation reduces corner separation loss obviously,positions from separation point to blade leading edge does not affect the final results. Combination actuation reduces loss of both suction surface and hub boundary layers to result in corner separation disappear and independence with attack angle.

corner separation;flow control;plasma actuation;phenomenological model;numerical simulation

梁斐杰（1987），男，工程師，在讀碩士研究生，研究方向?yàn)榈入x子體控制壓氣機(jī)角區(qū)分離流動(dòng)的數(shù)值模擬。

國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目（51136003）資助

2013-06-06