基于虛擬彈性體的快速動(dòng)網(wǎng)格方法

仲繼澤,徐自力,陶磊

(1.西安交通大學(xué)機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室,710049,西安;2.西安交通大學(xué)航天航空學(xué)院,710049, 西安; 3.國家電網(wǎng)浙江北侖第一發(fā)電有限責(zé)任公司,315800,浙江寧波)

?

基于虛擬彈性體的快速動(dòng)網(wǎng)格方法

仲繼澤1,2,徐自力1,2,陶磊3

(1.西安交通大學(xué)機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室,710049,西安;2.西安交通大學(xué)航天航空學(xué)院,710049, 西安; 3.國家電網(wǎng)浙江北侖第一發(fā)電有限責(zé)任公司,315800,浙江寧波)

為減少流固耦合的計(jì)算時(shí)間,在現(xiàn)有彈性體方法的基礎(chǔ)上,發(fā)展了一種快速動(dòng)網(wǎng)格方法。首先根據(jù)彈性體方法的基本假設(shè),將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,然后將結(jié)構(gòu)與該虛擬彈性體視為一個(gè)整體系統(tǒng),并計(jì)算其固有振動(dòng)的振型及頻率,最后以結(jié)構(gòu)受到的流體作用力為激勵(lì),通過振型疊加法計(jì)算結(jié)構(gòu)網(wǎng)格及流場網(wǎng)格節(jié)點(diǎn)的位移?？紤]到實(shí)際結(jié)構(gòu)的流固耦合振動(dòng)多為低階模態(tài)的振動(dòng),在流固耦合計(jì)算中可以通過低階模態(tài)的疊加計(jì)算流場網(wǎng)格節(jié)點(diǎn)的位移,從而達(dá)到快速更新流場網(wǎng)格的目的。采用該快速動(dòng)網(wǎng)格算法,對某彈性梁顫振問題進(jìn)行了流固耦合分析,計(jì)算結(jié)果與已有文獻(xiàn)的結(jié)果吻合很好,說明了該算法的正確性。與現(xiàn)有的彈性體方法相比,該算法使流固耦合計(jì)算時(shí)間減少了65.5%。對Wing 445.6模型的顫振問題進(jìn)行了分析,得到顫振邊界與實(shí)驗(yàn)值吻合良好,且與現(xiàn)有彈性體法相比,可以減少計(jì)算時(shí)間54.8%。

流固耦合;彈性體方法;快速動(dòng)網(wǎng)格方法

在流固耦合計(jì)算中,為考慮結(jié)構(gòu)振動(dòng)對流場的影響,通常需要采用動(dòng)網(wǎng)格方法更新流場網(wǎng)格[1]。目前所發(fā)展的網(wǎng)格變形方法主要有彈簧法[2]、彈性體方法[3]、溫度體方法[4]及徑向基函數(shù)方法[5]。與其他動(dòng)網(wǎng)格方法相比,彈性體方法是一種基于物理模型的網(wǎng)格變形方法,把流場網(wǎng)格變形作為一個(gè)彈性體變形問題進(jìn)行求解,更能貼近結(jié)構(gòu)彈性變形的物理實(shí)際。文獻(xiàn)[6]率先提出了彈性體方法,將流場網(wǎng)格所圍繞的空間區(qū)域視為虛擬彈性體,流場邊界變形作為彈性體的位移載荷,通過求解靜力平衡方程計(jì)算網(wǎng)格節(jié)點(diǎn)的位移。

在早期的研究中,文獻(xiàn)[6-8]都假定虛擬彈性體的彈性模量是各向同性的,采用這種彈性體方法更新流場網(wǎng)格時(shí),邊界層網(wǎng)格容易發(fā)生畸變,網(wǎng)格變形能力受到限制。文獻(xiàn)[9]引入了剛度系數(shù),將虛擬彈性體剛度與流場網(wǎng)格大小關(guān)聯(lián),通過增加小網(wǎng)格附近虛擬彈性體的剛度,使得網(wǎng)格變形能力和變形網(wǎng)格的質(zhì)量得到顯著提高。文獻(xiàn)[10]提出了分層彈性體方法,采用分層的思想,逐層決定虛擬彈性體的剛度,提高了網(wǎng)格變形能力。但是,在流固耦合計(jì)算中,每一時(shí)間步,都需要更新流場網(wǎng)格,且流場網(wǎng)格的數(shù)量通常為100萬甚至是1 000萬,因此流場網(wǎng)格變形的靜力平衡方程規(guī)模大,反復(fù)求解會(huì)占用大量的計(jì)算時(shí)間。

本文在彈性體方法的基礎(chǔ)上,發(fā)展了一種快速動(dòng)網(wǎng)格方法。將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,推導(dǎo)了結(jié)構(gòu)與虛擬彈性體同步振動(dòng)的動(dòng)力學(xué)方程。在流固耦合計(jì)算的每一時(shí)間步,通過低階模態(tài)的疊加,計(jì)算流場網(wǎng)格變形的節(jié)點(diǎn)位移,并更新流場網(wǎng)格。采用本文的動(dòng)網(wǎng)格算法,對彈性梁顫振算例進(jìn)行了流固耦合計(jì)算,計(jì)算結(jié)果與已有文獻(xiàn)的結(jié)果吻合很好,說明了本文算法的正確性。

1 快速動(dòng)網(wǎng)格方法

結(jié)構(gòu)的網(wǎng)格節(jié)點(diǎn)可以分成流固耦合面節(jié)點(diǎn)和非流固耦合面節(jié)點(diǎn)兩個(gè)部分,據(jù)此將結(jié)構(gòu)的振動(dòng)控制方程寫成分塊矩陣的形式

(1)

式中:Ms、Cs、Ks分別是結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣;Xs_in、Xs_n_in分別是結(jié)構(gòu)流固耦合面節(jié)點(diǎn)位移和非流固耦合面節(jié)點(diǎn)位移;Fin是結(jié)構(gòu)流固耦合面受到的流體作用力。

結(jié)構(gòu)振動(dòng)會(huì)改變流場的邊界,需要采用動(dòng)網(wǎng)格方法對流場網(wǎng)格節(jié)點(diǎn)坐標(biāo)進(jìn)行更新。在彈性體方法中,流場網(wǎng)格所圍繞的空間區(qū)域被視為虛擬彈性體,首先設(shè)定虛擬彈性體的彈性模量Ep,然后采用有限元方法計(jì)算得到的整體剛度矩陣Kp。在虛擬彈性體變形時(shí),流場網(wǎng)格節(jié)點(diǎn)位移Xp滿足靜力平衡方程[10]

(2)

式中:Xp_in、Xp_n_in分別是流場網(wǎng)格流固耦合面節(jié)點(diǎn)位移和非流固耦合面節(jié)點(diǎn)位移。對結(jié)構(gòu)的網(wǎng)格流固耦合面節(jié)點(diǎn)的位移Xs_in進(jìn)行插值,可以計(jì)算出流場網(wǎng)格流固耦合面節(jié)點(diǎn)的位移

Xp_in=NXs_in

(3)

式中:N是插值系數(shù)矩陣。流場網(wǎng)格變形的靜力平衡方程可以改寫為

(4)

將結(jié)構(gòu)與流場網(wǎng)格域的虛擬彈性體視為一個(gè)整體系統(tǒng),通過疊加結(jié)構(gòu)振動(dòng)控制方程和虛擬彈性體的靜力平衡方程,得到該整體系統(tǒng)的振動(dòng)方程

(5)

將該動(dòng)力學(xué)系統(tǒng)的固有頻率和正則振型分別記為ω和Ψ。式(5)經(jīng)坐標(biāo)變換后,得到正則坐標(biāo)ξ下的動(dòng)力學(xué)方程

(6)

式中:ζ是阻尼比;Q是正則坐標(biāo)下的流體作用力。

(7)

(8)

(9)

將式(6)寫成分量的形式如下

(10)

模態(tài)位移可以采用Wilson-θ方法[11]計(jì)算得到。實(shí)際結(jié)構(gòu)的流固耦合振動(dòng)多與低階模態(tài)相關(guān)[12],因此通過低階模態(tài)的疊加可以快速計(jì)算出結(jié)構(gòu)及流場網(wǎng)格節(jié)點(diǎn)的位移

(11)

式中:n是參與流固耦合振動(dòng)的模態(tài)的數(shù)量。

將固有頻率矩陣寫成分塊矩陣的形式

(12)

式中:ωs是結(jié)構(gòu)振動(dòng)主導(dǎo)的模態(tài)的固有頻率;ωp是虛擬彈性體振動(dòng)主導(dǎo)的模態(tài)的固有頻率。

固有模態(tài)是指振動(dòng)系統(tǒng)自由振動(dòng)時(shí)的固有頻率及振型。本文將結(jié)構(gòu)與流場網(wǎng)格域的虛擬彈性體作為一個(gè)整體系統(tǒng)。在整體系統(tǒng)的固有模態(tài)中,如果結(jié)構(gòu)部分基本保持了自身的固有模態(tài),將整體系統(tǒng)的固有模態(tài)稱為結(jié)構(gòu)振動(dòng)主導(dǎo)的模態(tài)。將振型矩陣寫成分塊矩陣的形式

(13)

式中:Ψss是整體系統(tǒng)振型對應(yīng)的結(jié)構(gòu)部分的振型;Ψpp是整體系統(tǒng)振型對應(yīng)的虛擬彈性體部分的振型。

根據(jù)固有頻率與剛度矩陣之間的關(guān)系式(8),結(jié)構(gòu)振動(dòng)主導(dǎo)的模態(tài)的固有頻率可以采用下式表達(dá)

(14)

將結(jié)構(gòu)部分的振動(dòng)Ψss代入到結(jié)構(gòu)的振動(dòng)控制方程式(1)中,可以得到結(jié)構(gòu)固有頻率的表達(dá)式為

(15)

與結(jié)構(gòu)的剛度相比,虛擬彈性體的剛度必須是小量。即,虛擬彈性體的彈性模量必須滿足下式

(16)

式中:Es為結(jié)構(gòu)的彈性模量。

2 算例驗(yàn)證

采用文獻(xiàn)[13]的彈性梁顫振算例驗(yàn)證本文提出的動(dòng)網(wǎng)格高效算法。流場及彈性梁的具體尺寸見圖1。梁的密度為1.0×104kg/m3,彈性模量為1.4×106Pa,泊松比為0.4。流體的密度為1 000 kg/m3,動(dòng)力學(xué)黏度為1 Pa·s。流場的左端面為流場的入口,采用速度進(jìn)口邊界條件,流速方向與進(jìn)口端面垂直,大小沿y方向呈拋物線分布

(17)

流場的右端面為流場的出口,采用壓力出口邊界,設(shè)定出口的平均靜壓值為0。流場的上端面和下端面、彈性梁的流固耦合面及左端的圓形壁面均采用無滑移邊界。

圖1 彈性梁及附近流場區(qū)域

流場及彈性梁的網(wǎng)格見圖2。梁的網(wǎng)格共有4節(jié)點(diǎn)矩形單元153個(gè),節(jié)點(diǎn)208個(gè);流場共有4節(jié)點(diǎn)矩形單元19 627個(gè),節(jié)點(diǎn)20 034個(gè)。

圖2 彈性梁及流場網(wǎng)格

2.1 流場網(wǎng)格域虛擬彈性體的彈性模量

本文將虛擬彈性體和結(jié)構(gòu)作為一個(gè)整體系統(tǒng)進(jìn)行研究,結(jié)構(gòu)的固有振動(dòng)特性勢必受到虛擬彈性體的影響。為了保證計(jì)算結(jié)果的正確性,必須使結(jié)構(gòu)基本上保持原有的固有振動(dòng)特性。如果設(shè)置的虛擬彈性體彈性模量不恰當(dāng),會(huì)使彈性梁振動(dòng)主導(dǎo)的模態(tài)頻率及振型偏離其固有振動(dòng)頻率和振型,從而產(chǎn)生錯(cuò)誤的計(jì)算結(jié)果,因此必須對流場網(wǎng)格設(shè)定合適的彈性模量。本文以彈性模量比對數(shù)為參數(shù),分析虛擬彈性體對彈性梁振動(dòng)主導(dǎo)的模態(tài)頻率及振型的影響,以確定虛擬彈性體的彈性模量。彈性模量比對數(shù)定義如下

R=lg(Es/Ep)

(18)

頻率偏差定義如下

(19)

式中:fs-p是彈性梁振動(dòng)主導(dǎo)的模態(tài)的固有頻率;fs是彈性梁的固有頻率。

頻率偏差與彈性模量比對數(shù)的關(guān)系見圖3。從圖中可以看出:在虛擬彈性體的彈性模量一定的情況下,隨著模態(tài)階數(shù)的增加,頻率偏差會(huì)下降;隨著彈性模量比對數(shù)的增加,頻率偏差逐漸減小;當(dāng)彈性模量比對數(shù)大于等于6,虛擬彈性體對彈性梁振動(dòng)主導(dǎo)的模態(tài)頻率的影響基本上可以忽略。

圖3 頻率偏差與彈性模量比對數(shù)的關(guān)系

彈性模量比對數(shù)與彈性梁振動(dòng)主導(dǎo)的模態(tài)振型之間的關(guān)系見圖4,其中位移和彈性梁長度均為無量綱量。隨著彈性模量比對數(shù)的增大,虛擬彈性體對彈性梁振型的影響逐漸減小。當(dāng)彈性模量比對數(shù)大于等于5,虛擬彈性體對彈性梁振動(dòng)主導(dǎo)的模態(tài)振型的影響基本上可以忽略。

(a)第1階振型

(b)第2階振型

(c)第3階振型

本文的研究結(jié)果表明,當(dāng)虛擬彈性體的彈性模量小于等于結(jié)構(gòu)彈性模量的1.0×106倍時(shí),虛擬彈性體對彈性梁固有頻率及振型的影響均可忽略。如果虛擬彈性體彈性模量取值過小,會(huì)超出計(jì)算機(jī)的精度范圍,造成浮點(diǎn)溢出。因此,本文建議虛擬彈性體彈性模量的取值如下

(d)第4階振型圖4 模態(tài)振型與彈性模量比對數(shù)的關(guān)系

(20)

此時(shí),彈性梁及虛擬彈性體整體系統(tǒng)的各階固有模態(tài)的振型見圖5。

(a)第1階模態(tài)

(b)第2階模態(tài)

(c)第3階模態(tài)

(d)第4階模態(tài)圖5 不同模態(tài)下的流場網(wǎng)格變形

2.2 流固耦合響應(yīng)分析

本文以定常流場的結(jié)果作為瞬態(tài)流場的初值,采用該動(dòng)網(wǎng)格算法對彈性梁顫振問題進(jìn)行分析。設(shè)定時(shí)間步長為0.001 s,流場變量的相對殘差隨迭代步變化的曲線見圖6。流場變量的相對殘差減小,在進(jìn)行25步迭代之后,殘差基本保持不變。

圖6 流場變量相對殘差與迭代次數(shù)的關(guān)系

(a)第1階

(b)第2階

(c)第3階

(d)第4階圖7 彈性梁的模態(tài)位移時(shí)間曲線

本文所采用的流動(dòng)條件正是文獻(xiàn)[13]中彈性梁顫振的流動(dòng)條件。在此流動(dòng)條件下,本文計(jì)算得到的彈性梁振動(dòng)的模態(tài)位移與時(shí)間曲線見圖7,模態(tài)位移為質(zhì)量歸一化量。隨著時(shí)間的推進(jìn),第1階,第3階和第4階振動(dòng)的模態(tài)位移幅值逐漸減小,即不會(huì)發(fā)生顫振。第2階模態(tài)振動(dòng)的位移幅值不隨時(shí)間變化,即第2階模態(tài)的振動(dòng)處于顫振臨界點(diǎn)。該彈性梁的顫振為第2階彎曲顫振。

(d)T圖8 振動(dòng)彈性梁附近的流場

采用流固耦合的方法,計(jì)算得到彈性梁右端中點(diǎn)的位移和時(shí)間曲線見圖9,從中可以看出,本文的計(jì)算結(jié)果與文獻(xiàn)[13]的結(jié)果吻合很好,驗(yàn)證了本文動(dòng)網(wǎng)格算法的正確性。本文計(jì)算得到圖9所示的彈性梁顫振的位移時(shí)間曲線,總耗時(shí)為19.8

h

。在相同計(jì)算條件下,采用傳統(tǒng)彈性體網(wǎng)格變形方法

[10]

的計(jì)算時(shí)間為57.4

h

,而本文方法可以減少計(jì)算時(shí)間65.5%。

為進(jìn)一步驗(yàn)證本文算法,我們對Wing445.6模型[14]的顫振問題進(jìn)行了流固耦合分析,計(jì)算得到的顫振邊界與實(shí)驗(yàn)值吻合良好,計(jì)算結(jié)果見圖10,與已有彈性體方法相比,計(jì)算時(shí)間減少了54.8%。

(a)x方向位移

(b)y方向位移圖9 彈性梁右端中點(diǎn)的位移和時(shí)間的關(guān)系

圖10 Wing 445.6模型的顫振邊界

3 結(jié) 論

本文提出了一種基于虛擬彈性體的快速動(dòng)網(wǎng)格方法。將流場網(wǎng)格所包圍的空間區(qū)域視為虛擬彈性體,并與結(jié)構(gòu)一起,視為一個(gè)整體系統(tǒng),以結(jié)構(gòu)受到的流體作用力為激勵(lì),通過振型疊加法計(jì)算結(jié)構(gòu)網(wǎng)格及流場網(wǎng)格節(jié)點(diǎn)的位移?？紤]到實(shí)際結(jié)構(gòu)的流固耦合振動(dòng)多為低階模態(tài)的振動(dòng),在流固耦合計(jì)算中可以通過低階模態(tài)的疊加計(jì)算流場網(wǎng)格節(jié)點(diǎn)的位移,從而達(dá)到快速更新流場網(wǎng)格的目的。采用該方法對彈性梁顫振問題進(jìn)行了流固耦合分析,計(jì)算結(jié)果與文獻(xiàn)的結(jié)果吻合很好。與已有彈性體方法相比,本文算法可以減少合計(jì)算時(shí)間65.5%。研究發(fā)現(xiàn),流場渦旋的周期性形成、發(fā)展及脫落是彈性梁顫振的主要原因。對Wing445.6模型顫振問題進(jìn)行了分析,得到顫振邊界與實(shí)驗(yàn)值吻合良好,且與已有彈性體法相比,可以減少計(jì)算時(shí)間54.8%。

[1] 張偉偉, 高傳強(qiáng), 葉正寅. 氣動(dòng)彈性計(jì)算中網(wǎng)格變形方法研究進(jìn)展 [J]. 航空學(xué)報(bào), 2014, 35(2): 303-319. ZHANG Weiwei, GAO Chuanqiang, YE Zhengyin. Progress on mesh deformation in computational aeroelasticity [J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(2): 303-319.

[2] 史愛明, 楊青, 楊永年. 非結(jié)構(gòu)運(yùn)動(dòng)網(wǎng)格下的三維機(jī)翼顫振數(shù)值分析 [J]. 振動(dòng)與沖擊, 2006, 24(6): 27-28. SHI Aiming, YANG Qing, YANG Yongnian. Numerical flutter analysis of a 3-D wing using unstructured dynamic mesh Euler method [J]. Journal of Vibration and Shock, 2006, 24(6): 27-28.

[3] TEZDUYAR T E. Stabilized finite element formulations for incompressible flow computations [J]. Advances in Applied Mechanics, 1991, 28: 1-44.

[4] 陳炎, 曹樹良, 梁開洪, 等. 基于溫度體模型的動(dòng)網(wǎng)格生成方法及在流固耦合振動(dòng)中的應(yīng)用 [J]. 振動(dòng)與沖擊, 2010, 29(4): 1-5. CHEN Yan, CAO Shuliang, LIANG Kaihong, et al. A new dynamic grids based on temperature analogy and its application in vibration engineering with fluid-solid interaction [J]. Journal of Vibration and Shock, 2010, 29(4): 1-5.

[5] 謝亮, 徐敏, 張斌, 等. 基于徑向基函數(shù)的高效網(wǎng)格變形算法研究 [J]. 振動(dòng)與沖擊, 2013, 32(10): 141-145. XIE Liang, XU Min, ZHANG Bin, et al. Space points reduction in grid deforming method based on radial basis functions [J]. Journal of Vibration and Shock, 2013, 32(10): 141-145.

[6] TEZDUYAR T E, BEHR M, MITTAL S, et al. A new strategy for finite element computations involving moving boundaries and interfaces-the deforming-spatial-domain/space-time procedure: I The concept and the preliminary tests [J]. Computer Methods in Applied Mechanics and Engineering, 1992, 94(3): 339-351.

[7] TEZDUYAR T E, BEHR M, MITTAL S, et al. A new strategy for finite element computations involving moving boundaries and interfaces-the deforming-spatial-domain/space-time procedure: II Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders [J]. Computer Methods in Applied Mechanics and Engineering, 1992, 94(3): 353-371.

[8] BAR-YOSEPH P Z, MEREU S, CHIPPADA S, et al. Automatic monitoring of element shape quality in 2-D and 3-D computational mesh dynamics [J]. Computational Mechanics, 2001, 27(5): 378-395.

[9] STEIN K, TEZDUYAR T, BENNEY R. Mesh moving techniques for fluid-structure interactions with large displacements [J]. Journal of Applied Mechanics, 2003, 70(1): 58-63.

[10]HUO S H, WANG F S, YAN W Z, et al. Layered elastic solid method for the generation of unstructured dynamic mesh [J]. Finite Elements in Analysis and Design, 2010, 46(10): 949-955.

[11]ZIENKIEWICZ O C, PAUL D K, HINTON E. Cavitation in fluid-structure response with particular reference to dams under earthquake loading [J]. Earthquake Engineering & Structural Dynamics, 1983, 11(4): 463-481.

[12]MARSHALL J G, IMREGUN M. A review of aeroelasticity methods with emphasis on turbomachinery applications [J]. Journal of Fluids and Structures, 1996, 10(3): 237-267.

[13]TUREK S, HRON J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow [M]. Berlin, Germany: Springer, 2006: 371-385.

[14]YATES E C, Jr. AGARD standard aeroelastic configurations for dynamic response I-wing 445.6 [R]. Neuilly Sur Seine, France: Advisory Group for Aerospace Research and Development Neuilly-Sur-Seine, 1988: 1-74.

(編輯 杜秀杰 趙煒)

An Efficient Dynamic Mesh Method Based on Pseudo Elastic Solid

ZHONG Jize1,2,XU Zili1,2,TAO Lei3

(1. State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University,Xi’an 710049, China; 2. School of Aerospace, Xi’an Jiaotong University, Xi’an 710049, China;3. State Grid Zhejiang Beilun No.1 Power Generation Co. Ltd., Ningbo, Zhejiang 315800, China)

To reduce the computing time of fluid structure coupling, an efficient dynamic mesh method based on the pre-existing elastic solid method is developed. The flow mesh domain is assumed to be a pseudo elastic solid according to the basic hypotheses of the elastic solid method, then the structure and the pseudo elastic solid are considered together as one holistic system. Subsequently the natural frequencies and vibration modes are calculated for the system and the flow force acted on the structure is considered as the excitation of the holistic system. The nodal displacements for the structure and the flow mesh are computed by mode superposition. In fact, the actual fluid structure coupled vibration for structures often appears associated with low order modes, the nodal displacements of the flow mesh can be calculated by modal superposition of the first few of low order modes, thus the flow mesh can be updated efficiently. A beam flutter problem is discussed with the present dynamic mesh method. The results coincide well with the data reported in the existing reference verifying the validation of the present method. The computing time is reduced by 65.5% compared with the pre-existing elastic solid method. The flutter of wing 445.6 is also analyzed and the calculated flutter boundary agrees with the experimental data, and the computing time is reduced by 54.8%.

fluid structure interaction; elastic solid method; efficient dynamic mesh method

2016-03-22。

仲繼澤(1988—),男,博士生;徐自力(通信作者),男,教授,博士生導(dǎo)師。

國家自然科學(xué)基金資助項(xiàng)目(51275385)。

時(shí)間:2016-07-14

http:∥www.cnki.net/kcms/detail/61.1069.T.20160714.1719.002.html

10.7652/xjtuxb201610020

O323;O351.2

A

0253-987X(2016)10-0132-07