基于正交檢驗法的顫振分析模態(tài)跟蹤技術(shù)

蔣令聞，杭曉晨，費慶國

(1. 江蘇省工程力學(xué)分析重點實驗室，南京 210096; 2. 東南大學(xué) 工程力學(xué)系，南京 210096)

基于正交檢驗法的顫振分析模態(tài)跟蹤技術(shù)

蔣令聞1,2，杭曉晨1,2，費慶國1,2

(1. 江蘇省工程力學(xué)分析重點實驗室，南京 210096; 2. 東南大學(xué) 工程力學(xué)系，南京 210096)

模態(tài)跟蹤是顫振分析的關(guān)鍵問題之一。常見的基于模態(tài)置信度的跟蹤方法在顫振點附近常會出現(xiàn)錯誤導(dǎo)致跟蹤失敗。通過探究該方法的數(shù)學(xué)機理，闡明了基于模態(tài)置信度跟蹤方法存在的缺陷，進而引入左特征向量，研究了一種基于正交檢驗法的模態(tài)跟蹤技術(shù)。采用一懸臂直翼盒模型開展仿真研究，在狀態(tài)空間下進行顫振分析，基于正交檢驗法得到其根軌跡圖，并將模態(tài)跟蹤結(jié)果與MAC值法進行對比，證明了正交檢驗法具有更好的模態(tài)跟蹤效果。

顫振；狀態(tài)空間；特征向量；模態(tài)跟蹤；模態(tài)置信度

顫振是飛行器在飛行中受到來流激勵后產(chǎn)生的一種災(zāi)難性的破壞振蕩，需要通過理論分析和設(shè)計來避免。通常基于結(jié)構(gòu)的各階模態(tài)頻率和阻尼隨空速的變化趨勢(即v-f圖和v-g圖)來預(yù)測顫振臨界點。確保模態(tài)在空速變化過程中的前后對應(yīng)關(guān)系是顫振分析的基礎(chǔ)，否則可能會給氣彈穩(wěn)定性判斷、顫振模態(tài)判別以及顫振抑制設(shè)計帶來困難。模態(tài)跟蹤技術(shù)旨在保證結(jié)構(gòu)的各階模態(tài)在系統(tǒng)參數(shù)變化過程中始終保持正確的前后對應(yīng)關(guān)系，不發(fā)生所謂的“模態(tài)交叉”現(xiàn)象[1]。因此，發(fā)展穩(wěn)定高效的模態(tài)跟蹤技術(shù)來防止顫振分析中的“模態(tài)交叉”是十分必要的。

顫振分析從數(shù)學(xué)角度上看是一個空速變化時的特征值問題，而模態(tài)跟蹤本質(zhì)上是對不同空速下的特征值進行匹配，從而使前后空速下得到的特征值保持正確的對應(yīng)關(guān)系。早期，分析員通常采用人工排序的方法，根據(jù)經(jīng)驗判斷前后空速下特征值的對應(yīng)關(guān)系，從而得到相對光滑合理的特征值曲線。這樣的做法非常繁瑣，而且結(jié)果的準(zhǔn)確性也難以保證。為此，國內(nèi)外研究人員發(fā)展了幾種模態(tài)跟蹤技術(shù)，包括：① 利用前后特征值之間的連續(xù)性對特征值進行排序[2]，這種方法在步長較長或模態(tài)較密集的情況下效果不好；② 基于特征值攝動理論的特征值預(yù)測方法[3-5]，即利用矩陣的左右特征向量信息計算特征值相對于空速的導(dǎo)數(shù)，從而給出下一空速的特征值預(yù)測值，以實現(xiàn)模態(tài)跟蹤。研究表明該方法的精度較高，但具有較大的計算量；③ 利用系統(tǒng)矩陣右特征向量來計算模態(tài)置信度(Modal Assurance Criterion，MAC值)[6-7]，從而實現(xiàn)模態(tài)跟蹤，此方法雖然簡單實用，但經(jīng)常在顫振點附近出現(xiàn)跟蹤錯誤，即使縮小步長也無法解決。

本文首先探究了傳統(tǒng)的MAC值法在應(yīng)用時出現(xiàn)失效的原因，從數(shù)學(xué)機理上明確了其存在的缺陷，并研究了一種基于正交檢驗法的模態(tài)跟蹤技術(shù)。采用懸臂直翼盒模型進行仿真研究，顫振分析結(jié)果表明，正交檢驗法相比于傳統(tǒng)的MAC值法更為有效。

1 理論方法

1.1氣彈運動方程與基于MAC的模態(tài)跟蹤技術(shù)

在模態(tài)坐標(biāo)下，線性氣動彈性系統(tǒng)的顫振運動方程可以寫為[8]

(1)

式中：M、D、K均為模態(tài)坐標(biāo)q下的廣義質(zhì)量矩陣、廣義阻尼矩陣和廣義剛度矩陣；qd=ρV2/2為空氣動壓；Q為廣義非定常氣動力系數(shù)矩陣，往往是頻域下的離散形式。

在氣彈伺服控制領(lǐng)域，通常在狀態(tài)空間模型下進行顫振分析。因此需要獲得拉氏域下的非定常氣動力矩陣。工程中常用有理多項式擬合的方法將頻域氣動力延拓至拉氏域。采用Roger法擬合多項式[9]代入式(1)并重寫為狀態(tài)空間的形式

(2)

式(2)可簡寫為

(3)

式(3)的特征值實部表征系統(tǒng)振動衰減率，虛部表征系統(tǒng)振動頻率。畫出特征值隨空速變化的根軌跡圖即可分析各空速下的穩(wěn)定性。由于x中包含了增量狀態(tài)空間坐標(biāo)qam，即系統(tǒng)矩陣發(fā)生了“擴階”，使式(3)中多出許多特征值，稱其為增根。模態(tài)跟蹤時，需要將它們與結(jié)構(gòu)模態(tài)區(qū)分開來。

MAC值法是工程中廣泛使用的一種模態(tài)跟蹤技術(shù)，其本質(zhì)是利用歸一化后模態(tài)振型的正交性來判斷兩模態(tài)是否一致。為了減少分析自由度，通常將問題轉(zhuǎn)化到模態(tài)坐標(biāo)系下考慮，此時其特征向量之間不僅關(guān)于質(zhì)量和剛度矩陣具有加權(quán)正交性，而且具有自正交性。故其特征向量間的MAC值矩陣的計算方法簡化為(本節(jié)所提到的特征向量均為右特征向量)

(4)

式中，ψi-1和ψi為前后空速下的特征向量以列向量形式組成的兩個矩陣。

將MAC值法應(yīng)用于狀態(tài)空間模型下的顫振分析，其主要思路為：根據(jù)式(4)計算前后空速的MAC值矩陣Si-1,i(由于增根的干擾，該矩陣一般不為方陣)。找到Si-1,i中每一行的最大值，認為其相關(guān)的兩階模態(tài)是正確對應(yīng)的。

MAC值法在模態(tài)跟蹤的前期往往效果較好，但在接近顫振速度時還是容易出現(xiàn)“模態(tài)交叉”現(xiàn)象。這是由于顫振分析的研究對象是氣動彈性系統(tǒng)，氣動力與結(jié)構(gòu)的運動相互耦合。相應(yīng)地，系統(tǒng)矩陣中就包含氣動力的成分，如氣動質(zhì)量、氣動剛度、氣動阻尼等。這些氣動力項是非對稱的[10]，導(dǎo)致原本對稱的質(zhì)量、剛度、阻尼矩陣不再對稱，從而使系統(tǒng)矩陣失去對稱性。根據(jù)特征值問題的性質(zhì)，這意味著其特征向量系同時失去了關(guān)于質(zhì)量、剛度矩陣的加權(quán)正交性和自正交性[11]。因此，基于式(4)所進行的模態(tài)跟蹤方法可能會失效，這點將在算例中予以驗證。

1.2正交檢驗法

針對MAC值法的缺陷，引入系統(tǒng)矩陣的左特征向量系φ來解決系統(tǒng)右特征向量系不正交的問題。調(diào)整矩陣S的算法如下

(5)

式中,φi-1為前一空速下系統(tǒng)矩陣的左特征向量矩陣，同樣以列向量的形式排列。新的系數(shù)矩陣Λi-1,i即為前后空速左右特征向量矩陣的乘積。可以證明，無論系統(tǒng)矩陣是否對稱，Λi,i在數(shù)學(xué)上是一個嚴(yán)格的對角矩陣，即系統(tǒng)矩陣的左右特征向量系是正交的。為了說明這一點，下面給出簡單的推導(dǎo)。

任意可對角化的方陣X可以被分解為

X=QEQ-1

(6)

式中,E為對角矩陣，若將式(6)兩邊同時右乘Q，化為

XQ=QE

(7)

這是一個標(biāo)準(zhǔn)右特征值問題的形式，E對角線上的元素為右特征值，Q的對應(yīng)列為相應(yīng)的右特征向量。同理，若將式(6)兩邊同時左乘Q-1，化為

Q-1X=EQ-1

(8)

這是一個標(biāo)準(zhǔn)左特征值問題的形式，它與右特征值問題具有相同的特征值，Q-1的每一行為其左特征向量?？梢?，矩陣X的左右特征向量矩陣可以化為互逆的形式。同時也說明矩陣左右特征向量系相互正交。特別的，當(dāng)X為一實對稱矩陣時，有

QTX=EQT

(9)

此時，X的左右特征向量恰好相等，于是右特征向量系自身即可滿足正交性。由于模態(tài)坐標(biāo)系下，M和K均為對角矩陣，轉(zhuǎn)化為標(biāo)準(zhǔn)特征值問題后所得的系統(tǒng)矩陣M-1K或K-1M均為對稱矩陣，故一般問題均可認為其右特征向量系是正交的。

然而，氣彈問題的系統(tǒng)矩陣由于氣動力項的影響，是非對稱的。因此，由于MAC值法只使用前后風(fēng)速下的右特征向量進行檢驗，其MAC值矩陣會喪失正交性。采用引入左特征向量的正交檢驗法后，新的系數(shù)矩陣Λi-1,i相較于Si-1,i具有更明顯的對角占優(yōu)性，特別是在顫振點附近，優(yōu)勢更為明顯。

圖1 正交檢驗法的流程圖Fig.1 Flowchart of the orthogonality check technique

此外，提出一種高效的系數(shù)矩陣計算方法。交換式(5)中左右特征向量矩陣的位置是完全可行的，即

(10)

這樣，只需間隔計算各空速下的左右特征向量即可完成模態(tài)跟蹤，對比原來的式(5)，減少了一半的計算量。正交檢驗法的具體流程與MAC值法相似，其流程圖如圖1所示。

2 仿真研究

2.1計算模型

本文選用一懸臂等截面直翼盒作為仿真模型，其幾何模型示意圖如圖2所示。該模型被很多研究者作為驗證模型使用，其詳細參數(shù)見文獻[12]。

翼盒的有限元結(jié)構(gòu)建模及其氣動力網(wǎng)格劃分如圖3所示。整個模型分為3個翼段，每個翼段包含兩塊蒙皮、兩塊梁腹板和一段翼肋，均設(shè)為殼單元；每個翼段的4段梁緣設(shè)為桿單元。模型材料設(shè)為鋁。利用Guyan減縮消除面內(nèi)位移，以此提高顫振分析的收斂性。前6階固有頻率與文獻值的對比如表1所示。模態(tài)頻率和振型均與文獻對應(yīng)的很好。

圖2 翼盒幾何模型(去掉上蒙皮)Fig.2 Geometric model of the wing box (upper skin removed)

2.2計算結(jié)果

在Nastran中使用偶極子格網(wǎng)法計算不同減縮頻率下的非定常氣動力，并采用表面樣條插值作用于結(jié)構(gòu)。不考慮結(jié)構(gòu)阻尼，從Nastran中提取計算所得的M、K和Q矩陣。采用Roger法得到擬合的氣動力有理多項式，從而將氣動彈性方程轉(zhuǎn)化為形如式(2)的狀態(tài)空間形式。通過對式(3)中系統(tǒng)矩陣H進行特征值求解可以畫出其根軌跡圖，從而判定系統(tǒng)穩(wěn)定性。一般認為當(dāng)某一階模態(tài)特征值實部發(fā)生“穿零”，即由負變正時，系統(tǒng)失去穩(wěn)定性。

圖3 翼盒有限元模型及氣動力模型Fig.3 FEM structural & aerodynamic model of the wingbox

表1 有限元模型前六階固有頻率與文獻值的對比Tab.1 The first six natural frequencies of FEM model and the reference values

如果不使用任何模態(tài)跟蹤方法，只保留虛部為正的特征值，并按虛部從小到大的順序取前六個特征值作為結(jié)構(gòu)模態(tài)，取步長為20 m/s，空速范圍為0～450 m/s，分別畫出其實部和虛部的根軌跡圖如圖4、圖5所示?？梢钥闯?，增根的出現(xiàn)產(chǎn)生了嚴(yán)重的干擾，導(dǎo)致V0之后跟蹤到的特征值均為增根，而不對應(yīng)真實的結(jié)構(gòu)模態(tài)。p-k法或者V-g法此類頻域方法，由于不受系統(tǒng)“擴階”的影響，計算得到的特征值總是對應(yīng)著某階結(jié)構(gòu)模態(tài)，也就是說，模態(tài)跟蹤只要將所有點正確的連接起來即可。而在狀態(tài)空間模型下，增根的數(shù)量往往是結(jié)構(gòu)模態(tài)的特征值的數(shù)倍，如果全部畫出這些點，肉眼根本無法分辨。因此，可靠高效的模態(tài)跟蹤技術(shù)對根軌跡法來說是不可或缺的。

在相同的輸入?yún)?shù)下，分別用MAC值法和正交檢驗法作出該系統(tǒng)的根軌跡圖，如圖6～圖9所示。MAC值法和正交檢驗法在2～6階模態(tài)的跟蹤上是一致的，由于2階扭轉(zhuǎn)模態(tài)實部“穿零”發(fā)生顫振，顫振速度大約為273 m/s，與文獻[13]中的結(jié)果267 m/s基本一致。而在第1階模態(tài)的跟蹤上，MAC值法出現(xiàn)了錯誤，在280 m/s處開始與第2階重合。事實上，第1階模態(tài)會在顫振速度之后頻率降為0，即出現(xiàn)發(fā)散現(xiàn)象，雖然發(fā)散速度在這里不是臨界速度。文獻[13-14]均給出了該發(fā)散速度的值，分別為287 m/s和292 m/s。顯然，MAC值法無法給出發(fā)散速度。從圖7、圖8可知，正交檢驗法成功地在顫振速度之后對1階模態(tài)進行跟蹤，發(fā)散速度為301 m/s。同時，與文獻[14]中的根軌跡圖進行對比，驗證了正交檢驗法在該例中的準(zhǔn)確性。

圖4 特征值虛部隨空速變化曲線(未使用模態(tài)跟蹤方法， ΔV=20 m/s)Fig.4 Velocity vs frequencies plot (without mode tracking, ΔV=20 m/s)

圖5 特征值實部隨空速變化曲線 (未使用模態(tài)跟蹤方法， ΔV=20 m/s)Fig.5 Velocity vs decay rate plot (without mode tracking, ΔV=20 m/s)

圖6 特征值虛部隨空速變化曲線(MAC值法，ΔV=20 m/s)Fig.6 Velocity vs frequencies plot (MAC method, ΔV=20 m/s)

圖7 特征值實部隨空速變化曲線 (MAC值法，ΔV=20 m/s)Fig.7 Velocity vs decay rate plot (MAC method, ΔV=20 m/s)

圖8 特征值虛部隨空速變化曲線(正交檢驗法，ΔV=20 m/s)Fig.8 Velocity vs frequencies plot (orthogonality check, ΔV=20 m/s)

圖9 特征值實部隨空速變化曲線 (正交檢驗法，ΔV=20 m/s)Fig.9 Velocity vs decay rate plot (orthogonality check, ΔV=20 m/s)

2.3對比與分析

雖然MAC值法在該例中出現(xiàn)的跟蹤錯誤并沒有影響臨界速度的判斷，但很難保證它在其他工程應(yīng)用，特別是在氣彈優(yōu)化問題中不出錯?？梢园l(fā)現(xiàn)，這種錯誤往往出現(xiàn)在臨近顫振點附近，且在發(fā)生模態(tài)重合現(xiàn)象的兩支或數(shù)支模態(tài)中。在本算例中，是由一階彎曲和二階扭轉(zhuǎn)耦合產(chǎn)生顫振，提取各空速下根據(jù)式(4)計算得到的1階、2階模態(tài)的MAC值Si,i(1,2)，畫出其隨空速的變化曲線，如圖10?？梢钥吹?，MAC值隨著空速逐漸變大，在顫振點附近達到了最大值，此后又逐漸回落。正如1.1節(jié)描述的那樣，氣動力項使各階模態(tài)的特征向量失去了正交性。對于兩階耦合模態(tài)而言，越靠近顫振點，其模態(tài)相似程度越高，在顫振點附近幾乎變成了同一階模態(tài)，因此其特征向量之間的夾角很小，MAC值接近于1。考慮到空速變化的影響，Si-1,i(1,2)的值在顫振點附近超過了Si-1,i(1,1)，使Si-1,i矩陣喪失了對角占優(yōu)性，導(dǎo)致了跟蹤失敗。

圖10 前兩階模態(tài)MAC值隨空速變化曲線Fig.10 Velocity vs the MAC value of the first two modes plot

為了更直觀的體現(xiàn)兩種方法在模態(tài)跟蹤時的表現(xiàn)，分別提取顫振點附近空速為280 m/s處，根據(jù)式(4)和式(5)計算得到的前六階模態(tài)之間的系數(shù)矩陣。將每行最大值歸一化后根據(jù)其元素的絕對值畫出三維柱狀圖，如圖11、圖12所示。從圖11可以看出，前兩階模態(tài)之間的MAC值與對角線元素大小基本相等。實際上我們通過圖6、圖7可知，此時Si-1,i(1,2)已經(jīng)超過了Si-1,i(1,1)，即矩陣喪失了對角占優(yōu)性。而從圖12來看，正交檢驗法獲得的系數(shù)矩陣對角占優(yōu)性依然十分明顯?？梢?，正交檢驗法在顫振點附近的穩(wěn)定性明顯高于MAC值法。

圖11 系數(shù)矩陣三維柱狀圖(MAC值法，V=280 m/s)Fig.11 3D bar graph of the coefficient matrix (MAC method, V=280 m/s)

圖12 系數(shù)矩陣三維柱狀圖(正交檢驗法，V=280 m/s)Fig.12 3D bar graph of the coefficient matrix (orthogonality check, V=280 m/s)

一般模態(tài)跟蹤方法可以考慮減小步長ΔV來避免出錯。而本算例中，將步長改為15 m/s甚至10 m/s，MAC值法依然會在顫振點附近失效，情況與圖6、圖7類似?？梢?，該算法在顫振點附近的收斂性也不好。而正交檢驗法在ΔV增至30 m/s時，效果依然穩(wěn)定，如圖13、圖14。由1.2節(jié)可知，同一空速下系數(shù)矩陣是嚴(yán)格的對角陣，因此即使出現(xiàn)跟蹤失敗，正交檢驗法通過適當(dāng)減小步長即可修正錯誤。

圖13 特征值虛部隨空速變化曲線(正交檢驗法，ΔV=30 m/s)Fig.13 Velocity vs frequencies plot (orthogonality check, ΔV=30 m/s)

圖14 特征值實部隨空速變化曲線 (正交檢驗法，ΔV=30 m/s)Fig.14 Velocity vs decay rate plot (orthogonality check, ΔV=30 m/s)

3 結(jié) 論

本文針對傳統(tǒng)MAC值法應(yīng)用于顫振分析時經(jīng)常出錯的問題，從數(shù)學(xué)機理角度解釋了其跟蹤失敗的原因，并研究了一種基于正交檢驗法的模態(tài)跟蹤技術(shù)。通過在狀態(tài)空間模型下的顫振分析算例對兩種方法進行對比，得出以下結(jié)論：

(1) MAC值矩陣在氣彈問題中容易在顫振點附近喪失對角占優(yōu)性，故MAC值法在顫振分析模態(tài)跟蹤中可能出錯，且很難通過縮小步長來解決。

(2) 正交檢驗法使系數(shù)矩陣的對角占優(yōu)性大大增加，對步長的要求較低，且對步長具有嚴(yán)格的收斂性。

(3) 從計算效率上來看，正交檢驗法只需要在每個空速下求解一次特征值問題，與MAC值法相比，在獲得更高的模態(tài)追蹤精度的同時，并未降低運算效率。

[1] ELDRED M S, VENKAYYA V B, ANDERSON W J. Mode tracking issues in structural optimization[J]. AIAA Journal, 1995, 33(10): 1926-1933.

[2] 謝長川, 胡薇薇, 楊超. 顫振分析中的特征值排序問題[J]. 數(shù)學(xué)的實踐與認識, 2007, 37(18): 141-146.

XIE Changchuan, HU Weiwei, YANG Chao. Eigenvalue sorting problem in flutter analysis[J]. Mathematics in Practice and Theory, 2007, 37(18): 141-146.

[3] 吳志剛, 宋晨, 楊超. 顫振分析中的模態(tài)跟蹤技術(shù)[J]. 北京航空航天大學(xué)學(xué)報, 2010，36(2): 163-167.

WU Zhigang, SONG Chen, YANG Chao. Mode tracking techniques in flutter analysis [J]. Journal of Beijing University of Aeronautics and Astronautics, 2010，36(2): 163-167.

[4] ELDRED M S, VENKAYYA V B, ANDERSON W J. New mode tracking methods in aeroelastic analysis[J]. AIAA Journal, 1995, 33(7): 1292-1299.

[5] 任智毅, 金海波, 丁運亮. 基于氣動力分段表達的顫振求解方法與模態(tài)跟蹤技術(shù)的應(yīng)用[J]. 空氣動力學(xué)學(xué)報, 2014, 32(2): 246-251.

REN Zhiyi, JIN Haibo, DING Yunliang. Flutter analysis based on the piecewise express of aerodynamic forces and the application of mode tracking technology[J]. Acta Aerodynamica Sinica, 2014, 32(2): 246-251.

[6] DESFORGES M J, COOPER J E, WRIGHT J R. Mode tracking during flutter testing using the modal assurance criterion[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 1996, 210(1): 27-37.

[7] GARIFFO J M. Generalized reduced order modeling of aeroservoelastic systems[D]. Los Angeles: University of California, 2013: 35-39.

[8] 管德. 飛機氣動彈性力學(xué)手冊[M]. 北京:航空工業(yè)出版社, 1994: 126-142.

[9] 趙永輝. 氣動彈性力學(xué)與控制[M]. 北京：科學(xué)出版社, 2007: 310-315.

[10] WRIGHT J R, COOPER J E, 姚一龍. 飛機氣動彈性力學(xué)及載荷導(dǎo)論[M]. 上海：上海交通大學(xué)出版社, 2010: 161-163.

[11] AFOLABI D，PIDAPARTI R M V，YANG H T Y. Flutter prediction using an eigenvector orientation approach[J]. AIAA Journal, 1998, 36(1): 69-74.

[12] BHATIA K G, RUDISILL C S. Optimization of complex structures to satisfy flutter requirements[J]. AIAA Journal, 1971, 9(8): 1487-1491.

[13] STRIZ A G, VENKAYYA V B. Influence of structural and aerodynamic modeling on flutter analysis[J]. Journal of Aircraft, 1994, 31(5): 1205-1211.

[14] HASAN M. Multidisciplinary design and optimization of a composite wing box[D]. Ankara:Middle East Technical University, 2003: 158-180.

Modetrackingtechniqueinflutteranalysisbasedonthemethodoforthogonalitycheck

JIANG Lingwen1, 2, HANG Xiaochen1, 2, FEI Qingguo1, 2

(1. Jiangsu Key Laboratory of Engineering Mechanics, Nanjing 210096, China;2. Department of Engineering Mechanics, Southeast University, Nanjing 210096, China)

Mode tracking is one of the key problems in flutter analysis. The common method based on the modal assurance criterion frequently fails near the critical point in its application. The mathematical mechanism of the common method was studied, thus revealing its defects in nature. A method based on the orthogonality check was proposed, in which the left eigenvector was introduced. The method was applied to the flutter analysis of an unswept cantilevered wing box in the state space and the root locus of the system was obtained. The numerical results were compared with those of the MAC method and it is validated that the orthogonality check is more effective in mode tracking.

flutter; state space; eigenvector; mode tracking; modal assurance criterion (MAC)

V211

A

10.13465/j.cnki.jvs.2017.19.013

國家自然科學(xué)基金(51408389;51438002);江蘇省自然科學(xué)基金青年基金項目(BK20140281)

基金項目：國家自然科學(xué)基金(11572086)；江蘇省研究生培養(yǎng)創(chuàng)新工程項目(KYLX15_0092)

2016-03-25 修改稿收到日期：2016-08-15

蔣令聞 男，博士，碩士生,1990年生

費慶國 男，教授，博士生導(dǎo)師,1977年生