于鵬垚,任慧龍,李陳峰,劉亞沖
(1.大連海事大學(xué) 交通運(yùn)輸裝備與海洋工程學(xué)院,遼寧 大連116026;2.哈爾濱工程大學(xué) 船舶工程學(xué)院,哈爾濱 150001)
船首外飄砰擊設(shè)計(jì)載荷直接計(jì)算
于鵬垚1,任慧龍2,李陳峰2,劉亞沖2
(1.大連海事大學(xué) 交通運(yùn)輸裝備與海洋工程學(xué)院,遼寧 大連116026;2.哈爾濱工程大學(xué) 船舶工程學(xué)院,哈爾濱 150001)
文章結(jié)合三維線性勢流理論和砰擊速度的長期分析方法,求解出船體外飄位置的設(shè)計(jì)砰擊速度;以首垂線和靜水面交點(diǎn)處的設(shè)計(jì)砰擊速度為目標(biāo)值,給出了用于確定船首外飄砰擊設(shè)計(jì)載荷的等效設(shè)計(jì)波,進(jìn)而得到了設(shè)計(jì)狀態(tài)下船體外飄剖面與波面相對運(yùn)動關(guān)系;將船體剖面與波面間的相對運(yùn)動關(guān)系等效轉(zhuǎn)化為船體剖面與靜水面的相對運(yùn)動,利用顯式有限元方法實(shí)現(xiàn)了外飄剖面砰擊設(shè)計(jì)載荷的預(yù)報(bào)。針對直接計(jì)算方法中涉及的設(shè)計(jì)砰擊速度、砰擊壓力和砰擊壓力系數(shù),對比分析了文中結(jié)果和相應(yīng)的規(guī)范值或試驗(yàn)值,論證了文中船舶外飄砰擊壓力設(shè)計(jì)載荷直接計(jì)算方法的合理性。
外飄砰擊載荷;設(shè)計(jì)砰擊速度;砰擊壓力
船舶在高海況航行時(shí),由于船體與波浪之間劇烈的沖擊作用,具有外飄形式的船首結(jié)構(gòu)往往承受較大的砰擊載荷。近年來,由砰擊載荷導(dǎo)致船首局部結(jié)構(gòu)損傷的案例也時(shí)有耳聞,因此許多船級社規(guī)范都明確提出了考慮砰擊載荷下的船首結(jié)構(gòu)設(shè)計(jì)尺寸要求。然而可能由于許多細(xì)節(jié)方面的處理方法不一致,導(dǎo)致不同船級社對船首外飄結(jié)構(gòu)尺寸要求存在差別[1]。外飄砰擊壓力設(shè)計(jì)載荷的準(zhǔn)確預(yù)報(bào)是合理地設(shè)計(jì)首部抗砰擊結(jié)構(gòu)尺寸的前提,因此,十分有必要形成一套完整的外飄砰擊設(shè)計(jì)載荷確定方法,從而為船級社推出“抗砰擊結(jié)構(gòu)設(shè)計(jì)的共同規(guī)范”提供基礎(chǔ)。
Zhao[2-3],Sun等[4]采用非線性邊界元法針對不同結(jié)構(gòu)形式的二維剖面砰擊載荷進(jìn)行預(yù)報(bào),并與相應(yīng)的試驗(yàn)值進(jìn)行對比,驗(yàn)證數(shù)值算法的準(zhǔn)確性;Hermundstad[5]結(jié)合非線性切片法和非線性邊界元法預(yù)報(bào)了波浪中船舶砰擊載荷,為實(shí)現(xiàn)計(jì)算,對含球鼻首結(jié)構(gòu)的橫剖面處進(jìn)行了光順處理;駱寒冰[6]和Wang[7]采用LS-DYNA軟件對自由下落楔形體的入水砰擊載荷進(jìn)行仿真,并與試驗(yàn)結(jié)果吻合較好。Veen[8]將SPH(Smoothed Particle Hydrodynamics)方法與切片法結(jié)合,預(yù)報(bào)了迎浪工況下船首的外飄砰擊載荷。上述文獻(xiàn)主要針對某一工況下,船舶結(jié)構(gòu)所遭受的砰擊載荷的數(shù)值預(yù)報(bào)方法進(jìn)行研究,船舶一生中可能遭受各種各樣的波浪,如何選取一個(gè)合理的工況來確定船體結(jié)構(gòu)的砰擊載荷設(shè)計(jì)值仍是設(shè)計(jì)者十分關(guān)心的問題。王輝[9]利用非線性切片理論和Stavovy-Chang砰擊理論針對波浪中的船舶外飄砰擊載荷進(jìn)行預(yù)報(bào),并通過對短期海況的砰擊壓力值的統(tǒng)計(jì)分析,給出了局部砰擊壓力載荷的設(shè)計(jì)值;田喜民[10]采用三維Rankine源法求解船體在不規(guī)則波中的運(yùn)動,進(jìn)而結(jié)合砰擊壓力系數(shù),給出砰擊壓力峰值,同樣通過短期分析的方法給出了砰擊載荷的設(shè)計(jì)值。然而,在采用短期方法進(jìn)行設(shè)計(jì)值的預(yù)報(bào)時(shí),如何選取合理的短期海況仍有待商榷,而且文獻(xiàn)[9-10]的方法僅給出砰擊壓力峰值,未給出結(jié)構(gòu)動力響應(yīng)分析時(shí)所需的砰擊壓力隨時(shí)間的變化關(guān)系。
借鑒相關(guān)文獻(xiàn)的研究思路,本文形成了一套船首外飄砰擊設(shè)計(jì)載荷計(jì)算方法。利用海浪的長期統(tǒng)計(jì)資料對首垂線與靜水面交點(diǎn)處的垂向相對速度進(jìn)行統(tǒng)計(jì)分析,選取所需超越概率水平的下設(shè)計(jì)值來確定設(shè)計(jì)波;將船體在此設(shè)計(jì)波下所遭受的局部砰擊壓力載荷作為設(shè)計(jì)值,利用設(shè)計(jì)波下的船波運(yùn)動關(guān)系確定了船首隨時(shí)間的變化外飄砰擊載荷。具體的外飄砰擊設(shè)計(jì)載荷直接計(jì)算過程如圖1。
圖1 外飄砰擊設(shè)計(jì)載荷直接計(jì)算流程圖Fig.1 The flowchart for the direct calculation of bowflare slamming design loads
1.1 船波相對運(yùn)動關(guān)系
針對某一固定海況而言,頻域理論和時(shí)域理論均已被廣泛地應(yīng)用于船波在波浪中運(yùn)動與載荷的預(yù)報(bào),但在通過長期方法確定船舶結(jié)構(gòu)設(shè)計(jì)載荷時(shí),計(jì)算相對方便的頻域方法被采用得更為廣泛,本文采用三維線性頻域勢流理論進(jìn)行船波相對運(yùn)動的預(yù)報(bào)。當(dāng)重心處輸入的單位波幅簡諧波eiωet激勵下,船舶的運(yùn)動響應(yīng)也可表達(dá)為如下形式:
考慮到船舶運(yùn)動和波面的升高,船體上任意一點(diǎn)xb,yb,zb的垂向相對運(yùn)動響應(yīng)表達(dá)為:
式中:Hz(ωe)為單位波幅下垂向相對運(yùn)動;ωe為遭遇頻率;{η }為船舶六自由度響應(yīng);η3,η4,η5分別為單位波幅波浪作用下船舶升沉運(yùn)動、橫搖運(yùn)動和縱搖運(yùn)動;k為波數(shù);β為浪向角。
在實(shí)際的船舶砰擊過程中,船波相對運(yùn)動由垂向相對運(yùn)動和水平相對運(yùn)動組成,垂向相對運(yùn)動起主要貢獻(xiàn),水平相對運(yùn)動貢獻(xiàn)較小。對于沿船長方向型線變化不是十分劇烈的外飄區(qū)域,參考文獻(xiàn)[11]的處理方法,僅考慮垂向相對運(yùn)動對砰擊速度的貢獻(xiàn)量。
1.2 砰擊速度長短期分析
對于波浪譜密度函數(shù)為Sω(ωe)的短期海況,可以得到該海況下相對位移的方差和相對速度的方差分別為:
船體局部位置砰擊次數(shù)的計(jì)算在數(shù)學(xué)上屬于隨機(jī)過程中的過閾問題,船體表面距靜水面為zi的局部位置的單位時(shí)間以砰擊速度zv<=zv′發(fā)生砰擊次數(shù)Ns可以表達(dá)為[12]:
其中:zv(t)為船波垂向相對運(yùn)動速度,發(fā)生砰擊時(shí)。
當(dāng)zi=zv′=0時(shí),可得船舶在單位時(shí)間內(nèi)遭遇波浪的次數(shù)N0為:
考慮到航速,航向角和海情的變化,參考文獻(xiàn)[11]給定的航速與海情相互對應(yīng)的關(guān)系,可以得到砰擊速度zv<=zv′時(shí)概率分布函數(shù)F(zv′):
其中:pi(H1/3,Tz)為有義波高為H1/3、跨零周期為Tz的海況出現(xiàn)的概率;pj(β)為航向角β出現(xiàn)的概率;Ns(i,j),N0(i,j)分別為相應(yīng)短期海況下,砰擊速度zv<=zv′時(shí)的單位時(shí)間砰擊次數(shù)和船舶單位時(shí)間遭遇波浪的次數(shù)。
1.3 設(shè)計(jì)狀態(tài)下船波運(yùn)動關(guān)系
在確定船舶的設(shè)計(jì)載荷時(shí),通常規(guī)定船舶一生中遭遇的波浪載荷循環(huán)次數(shù)n=108,則計(jì)算時(shí)可以取概率水平F(zv′)=1/n=10-8下的砰擊速度值作為砰擊速度設(shè)計(jì)值。對于砰擊速度設(shè)計(jì)值,其表達(dá)的含義為在船舶一生遭遇的108個(gè)波浪中,會有一個(gè)較大的波浪使其局部的砰擊速度達(dá)到相應(yīng)的砰擊速度設(shè)計(jì)值,而這樣的一個(gè)波浪通常理解為用于結(jié)構(gòu)強(qiáng)度設(shè)計(jì)的設(shè)計(jì)波。在一個(gè)設(shè)計(jì)波下,船首外飄是一個(gè)區(qū)域,并不能夠同時(shí)達(dá)到采用長期方法得到的砰擊載荷設(shè)計(jì)值,通過后續(xù)計(jì)算可以看出,局部位置越靠近船首,越靠近水線時(shí),其砰擊速度越大,因而本文取首垂線與靜水面交點(diǎn)處的砰擊速度設(shè)計(jì)值為目標(biāo)值。選取單位波幅下不同浪向和不同頻率下該位置相對速度較大的浪向和頻率作為設(shè)計(jì)波的浪向和頻率,進(jìn)而利用砰擊速度設(shè)計(jì)值便可確定設(shè)計(jì)波的波幅。確定設(shè)計(jì)波后,便可得到設(shè)計(jì)狀態(tài)下船體外飄剖面與波面之間的相對運(yùn)動關(guān)系。
本文采用了顯示有限元分析程序LS-DYNA進(jìn)行砰擊載荷的數(shù)值預(yù)報(bào)。在具體的數(shù)值計(jì)算過程中,水和空氣采用多物質(zhì)的歐拉網(wǎng)格進(jìn)行模擬,剖面結(jié)構(gòu)采用拉格朗日網(wǎng)格進(jìn)行模擬,自由面的生成和重構(gòu)采用有限體積法。通過在歐拉網(wǎng)格(流體)和拉格朗日網(wǎng)格(結(jié)構(gòu))表面引入罰函數(shù)耦合算法,來實(shí)現(xiàn)流體與固體之間的耦合效應(yīng)的模擬。為驗(yàn)證本文采用方法在波浪中外飄砰擊載荷預(yù)報(bào)的實(shí)用性,將采用顯示有限元方法的預(yù)報(bào)結(jié)果與文獻(xiàn)[8]的結(jié)果進(jìn)行對比。其中,t=0時(shí)刻為剖面最低點(diǎn)觸及靜水面的時(shí)刻。通過圖4和圖5不同的方法壓力計(jì)算結(jié)果的比較可以看出,本文所采用的方法能夠滿足工程應(yīng)用的計(jì)算精度要求。
圖2 船首外飄剖面Fig.2 Bowflare section
圖3 剖面與波面間的相對速度Fig.3 Relative velocity between the section and the wave surface
圖5 不同方法的P2位置壓力值Fig.5 Location P2pressure of different methods
針對某大型集裝箱船,進(jìn)行外飄砰擊壓力設(shè)計(jì)載荷的直接計(jì)算。具體的船舶參數(shù)(表1)及計(jì)算分析結(jié)果如下。
3.1 船舶參數(shù)
表1 主尺度參數(shù)Tab.1 Principal dimensions
3.2 砰擊速度設(shè)計(jì)值
針對船首外飄區(qū)域的關(guān)注位置和首垂線與靜水面的交點(diǎn)進(jìn)行砰擊速度的長期分析。其中,海浪資料選取為國際船級社協(xié)會(IACS)推薦的北大西樣海浪長期統(tǒng)計(jì)資料;長峰波波幅選取為國際力學(xué)船舶結(jié)構(gòu)會議(ISSC)建議的雙參數(shù)譜;引入擴(kuò)散函數(shù),將長峰波譜轉(zhuǎn)化為短峰波幅,其中θ為組合波與主浪向的夾角。選取超越概率為10-8的砰擊速度值作為上述位置砰擊速度的設(shè)計(jì)值。通過首垂線與靜水面交點(diǎn)的相對運(yùn)動確定船首砰擊載荷設(shè)計(jì)波的浪向?yàn)橛耍ɡ俗匀活l率為0.5 rad/s,利用該位置超越概率水平下的設(shè)計(jì)值,確定該設(shè)計(jì)波的波幅為11.62 m。通過此設(shè)計(jì)波也可確定外飄關(guān)注位置的砰擊速度,后續(xù)稱之為設(shè)計(jì)波法。同樣針對相同的外飄關(guān)注位置,分別采用ABS指南[11]和LR[13]規(guī)范進(jìn)行垂向砰擊速度的計(jì)算。不同方法的計(jì)算結(jié)果對比見表2和圖7。
表2 砰擊壓力計(jì)算點(diǎn)位置Tab.2 Locations for calculating slamming pressure
圖6 計(jì)算點(diǎn)位置示意圖Fig.6 The schematic of calculation location s
圖7 不同方法確定的設(shè)計(jì)砰擊速度Fig.7 The design impact velocity of different methods
表3 不同方法確定的設(shè)計(jì)砰擊速度值(m/s)Tab.3 Design impact velocity of different methods
通過不同方法確定的設(shè)計(jì)砰擊速度的比較(表3,圖6-7)可以看出,長期值、設(shè)計(jì)波值、ABS指南和LR規(guī)范的計(jì)算結(jié)果依次減小,但總體變化趨勢一致??紤]到設(shè)計(jì)波是本文方法通過合理的長期分析確定,LR規(guī)范只是通過一些簡化公式給出設(shè)計(jì)值,ABS指南短期海況的選取合理性仍有待論證,而且設(shè)計(jì)波法更加方便確定砰擊載荷時(shí)外飄剖面運(yùn)動軌跡的確定,故建議采用設(shè)計(jì)波法進(jìn)行設(shè)計(jì)狀態(tài)下船體剖面與波面相對運(yùn)動關(guān)系的確定。
3.3 砰擊載荷預(yù)報(bào)
利用確定預(yù)報(bào)船首外飄砰擊載荷的設(shè)計(jì)波,預(yù)報(bào)得到船體第18站和第19站外飄剖面的入水位移和入水速度(圖8-9),其中t=0時(shí)刻分別取為各自剖面最低點(diǎn)觸及波面的時(shí)刻。
圖8 外飄剖面入水砰擊的位移Fig.8 Relative displacement between the section and the wave surface
圖9 外飄剖面入水砰擊的速度Fig.9 Relative velocity between the section and the wave surface
圖10 t=1.1 s時(shí)18站剖面的自由液面形狀Fig.10 Free surface shape of section 18 at t=1.1 s
圖11 t=1.1 s時(shí)19站剖面的自由液面形狀Fig.11 Free surface shape of section 19 at t=1.1 s
將船體外飄剖面進(jìn)入波面的過程等價(jià)為船體剖面進(jìn)入靜水面的過程,從而實(shí)現(xiàn)船體剖面入水砰擊載荷的預(yù)報(bào)。如圖10-13所示。
圖12 18站剖面的砰擊壓力Fig.12 Slamming pressure of section 18
圖13 19站剖面的砰擊壓力Fig.13 Slamming pressure of section 19
圖14 直接計(jì)算確定的砰擊壓力系數(shù)與規(guī)范值的比較Fig.14 Comparison between the direct calculation method and rules
3.4 砰擊壓力峰值系數(shù)
在進(jìn)行砰擊壓力載荷的預(yù)報(bào)時(shí),已往文獻(xiàn)[14]經(jīng)常利用砰擊速度和經(jīng)驗(yàn)公式給定砰擊壓力系數(shù)換算得到,本文利用預(yù)報(bào)得到的砰擊壓力峰值和峰值時(shí)刻對應(yīng)的沖擊速度換算得到相應(yīng)的砰擊壓力系數(shù)(表4),并與國家軍用標(biāo)準(zhǔn)(GJB)[15]和LR規(guī)范[13]經(jīng)驗(yàn)公式結(jié)果的比較見圖14??梢钥闯觯鼺1901和F1902兩個(gè)測點(diǎn)外,其他測點(diǎn)均在LR規(guī)范曲線附近,主要是由于規(guī)范曲線通常是由針對不同底升角的楔形體入水試驗(yàn)或理論計(jì)算得到,而楔形體在入水過程中,不存類似于19站剖面的,先發(fā)生流動分離再砰擊到結(jié)構(gòu)表面的情況。通過本文的計(jì)算可以看出,當(dāng)流體在球鼻首處發(fā)生流動分離后,再沖擊到結(jié)構(gòu)表面時(shí)往往產(chǎn)生更大的壓力,更應(yīng)該引起結(jié)構(gòu)設(shè)計(jì)的關(guān)注。
表4 直接計(jì)算確定的砰擊壓力系數(shù)(m/s)Tab.4 Slamming pressure coefficient of the direct calculation method
本文通過對與外飄砰擊壓力設(shè)計(jì)載荷相關(guān)的設(shè)計(jì)砰擊速度和砰擊壓力載荷預(yù)報(bào)兩方面的研究,得到如下結(jié)論:
(1)結(jié)合耐波性理論和長期分析方法,本文采用設(shè)計(jì)波法合理地預(yù)報(bào)了設(shè)計(jì)狀態(tài)下外飄結(jié)構(gòu)的沖擊速度及與其與波面的相對運(yùn)動關(guān)系。相比于短期預(yù)報(bào)方法和經(jīng)驗(yàn)公式預(yù)報(bào)方法,本文方法避免了短期設(shè)計(jì)海況的選取和經(jīng)驗(yàn)公式可能無法反映出船型特點(diǎn)等問題。
(2)利用顯式有限元方法,本文實(shí)現(xiàn)了外飄剖面砰擊設(shè)計(jì)載荷的預(yù)報(bào),并與相關(guān)文獻(xiàn)進(jìn)行了對比驗(yàn)證。相比于只能夠預(yù)報(bào)砰擊壓力峰值的設(shè)計(jì)載荷預(yù)報(bào)方法,本文方法可以得到入水過程中的砰擊壓力變化量,更方便后續(xù)進(jìn)一步與結(jié)構(gòu)動力分析方法相結(jié)合,實(shí)現(xiàn)砰擊強(qiáng)度的評估。
(3)針對受球鼻首影響的外飄剖面,本文方法能夠反映出球鼻首處發(fā)生流動分離后,水流飛濺到外飄區(qū)域的砰擊問題,而基于楔形體入水理論的經(jīng)驗(yàn)公式方法并不能很好地解決這一問題。通過本文的計(jì)算可以看出飛濺的水流沖擊到結(jié)構(gòu)表面時(shí)往往產(chǎn)生更大的壓力,更應(yīng)該引起結(jié)構(gòu)設(shè)計(jì)的關(guān)注。
[1]ISSC.Committee V.7.Impulsive pressure loading and response assessment[C].17th International Ship and Offshore Structures Congress,2012:31-33.
[2]Zhao R,Faltinsen O M.Water entry of two-dimensional bodies[J].Journal of Fluid Mechanics,1991,222:215-230.
[3]Zhao R,Faltinsen O M,Aarsnes J V.Water entry of arbitrary two-dimensional sections with and without flow separation [C].Proc.Twenty-first Symposium on Naval Hydrodynamics,1996:408-423.
[4]Sun Hui.A boundary element method applied to strongly nonlinear wave-body interaction problems[D].Norway:Norwegian University of Science and Technology,2007:45-65.
[5]Hermundstad O A,Moan T.Numerical and experimental analysis of bow flare slamming on a Ro-Ro vessel in regular oblique waves[J].J Mar Sci Technol,2005,10:105-122.
[6]駱寒冰,吳景健,王 珊,徐 慧.基于顯式有限元方法的二維楔形剛體入水砰擊載荷并行計(jì)算預(yù)報(bào)[J].船舶力學(xué), 2012,16(8):907-913. Luo Hanbing,Wu Jingjian,Wang Shan,Xu Hui.Parallel computing simulation of water entry of a 2D rigid wedge using an explicit finite element method[J].Journal of Ship Mechanics,2012,16(8):907-914.
[7]Wang Shan.Assessment of slam induced loads on two dimensional wedges and ship sections[D].Portugal:University of Lisbon,2011:27-40.
[8]Veen D,Gourlay T.A combined strip theory and Smoothed Particle hydrodynamics approach for estimating slamming loads on a ship in head seas[J].Ocean Engineering,2012,43:64-71.
[9]王 輝,朱云祥,顧學(xué)康,沈進(jìn)威.Ro-Ro船首門砰擊載荷設(shè)計(jì)值的確定[J].船舶力學(xué),2003,7(1):46-54. Wang Hui,Zhu Yunxiang,Gu Xuekang,Shen Jinwei.Determination of design slamming loads on bow doors for Ro-Ro ships[J].Journal of Ship Mechanics,2003,7(1):46-54.
[10]田喜民,鄒早建,王福花.大型船舶外飄砰擊壓力計(jì)算研究[J].中國造船,2014,55(1):1-10. Tian Ximin,Zou Zaojian,Wang Fuhua.Numerical study on flare impact for large ships[J].Shipbuilding of China,2014, 55(1):1-10.
[11]American Bureau of Shipping.Guide for slamming loads and strength assessment for vessels[S].2011.
[12]楊代盛,桑國光,李維揚(yáng),戴仰山.船舶強(qiáng)度的概率方法[M].哈爾濱:哈爾濱工程大學(xué)出版社,2007:100-119. Yang Daisheng,Sang Guoguang,Li Weiyang,Dai Yangshan.Probability method of ship strength[M].Harbin:Harbin Engineering University Press,2007:100-119.
[13]Rules and regulations for the classification of ships[S].London:Lloyd’s Register,2012.
[14]王 輝.船體結(jié)構(gòu)局部強(qiáng)度設(shè)計(jì)中的砰擊載荷確定方法[J].中國造船,2010,51(2):68-77. Wang Hui.Slamming load determination in local structure design of ships[J].Shipbuilding of China,2010,51(2):68-77.
[15]中華人民共和國國家軍用標(biāo)準(zhǔn) GJB64.1-85,艦船船體規(guī)范—水面艦艇[S].國防科學(xué)技術(shù)工業(yè)委員會批準(zhǔn),1985-08-31發(fā)布. National Militrary Standards of People’s Republic of China GJB/Z-1999,surface ship structual design calculation methods [S].Defense Technology Industry Council,1999.
Direct calculation method of bowflare design slamming loads
YU Peng-yao1,REN Hui-long2,LI Chen-feng2,LIU Ya-chong2
(1.Transportation Equipment and Ocean Engineering College,Dalian Maritime University,Dalian 116026,China; 2.College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China)
The three-dimensional linear potential flow theory and long-term analysis method of relative slamming velocity are combined to determine the design impact velocity of bowflare structures;the design impact velocity of the location where the forward perpendicular and the static water surface intersect is set as the target value to calculate the equivalent design wave amplitude and the relative motion between the bowflare section and the wave surface in the design state;the relative motion between the hull flare section and the wave surface is equivalently transformed as that between the hull section and the static water surface;slamming pressure loads are gained by an explicit finite element method.The results such as design impact velocity,slamming pressure and slamming pressure coefficient are compared with those of the rule or experiments,and through analyzing the difference,the direct calculation method of bowflare design slamming pressure loads shows reasonable.
bowflare slamming loads;design impact velocity;slamming pressure
U661.3
:Adoi:10.3969/j.issn.1007-7294.2016.05.007
1007-7294(2016)05-0566-08
2015-12-29
國家973基金資助項(xiàng)目(2011CB013703)
于鵬垚(1988-),男,博士,講師,E-mail:pengyao_yu@126.com;任慧龍(1965-),男,教授,博士生導(dǎo)師,通信作者。