梁崗 曹沛 唐炯 沈晴 李心爽

(1.上海海事大學(xué)物流工程學(xué)院, 上海 201306) (2.上海振華重工(集團(tuán))股份有限公司, 上海 200125)(3.美國(guó)康涅狄格大學(xué)機(jī)械工程系, 美國(guó))

基于移動(dòng)有限元法的裂紋梁振動(dòng)分析

梁崗1曹沛2?唐炯3沈晴2李心爽2

(1.上海海事大學(xué)物流工程學(xué)院, 上海 201306) (2.上海振華重工(集團(tuán))股份有限公司, 上海 200125)(3.美國(guó)康涅狄格大學(xué)機(jī)械工程系, 美國(guó))

采用移動(dòng)有限元法和局部柔度法對(duì)移動(dòng)質(zhì)量作用下含裂紋簡(jiǎn)支梁進(jìn)行了振動(dòng)計(jì)算分析.計(jì)算考慮了裂紋和移動(dòng)質(zhì)量的相對(duì)位置對(duì)梁固有頻率的影響,以及移動(dòng)質(zhì)量在不同位置、速度情況下對(duì)裂紋梁的動(dòng)力響應(yīng)的影響.結(jié)果分析表明,裂紋與移動(dòng)質(zhì)量的存在會(huì)使得梁的動(dòng)態(tài)位移有不同程度的增大,且隨著移動(dòng)質(zhì)量位置和裂紋位置的改變會(huì)使得梁的固有頻率變小.

移動(dòng)有限元, 簡(jiǎn)支梁, 張開(kāi)裂紋, 固有頻率, 動(dòng)力響應(yīng)

引言

移動(dòng)載荷對(duì)梁結(jié)構(gòu)動(dòng)力特性的影響在港口起重運(yùn)輸機(jī)械上是一個(gè)很重要的課題,如橋式起重機(jī)載重運(yùn)行和岸邊集裝箱起重機(jī)裝卸作業(yè)等均可簡(jiǎn)化為移動(dòng)質(zhì)量-梁耦合系統(tǒng).目前,工程實(shí)踐中常以梁的固有頻率代替移動(dòng)質(zhì)量-梁的耦合頻率,但是考慮到移動(dòng)質(zhì)量在運(yùn)動(dòng)過(guò)程中,質(zhì)量-梁組成的系統(tǒng)會(huì)產(chǎn)生耦合振動(dòng),對(duì)梁的固有頻率產(chǎn)生較大影響；另外,由于各種因素的影響,梁結(jié)構(gòu)中極易產(chǎn)生微觀裂紋.裂紋的存在不僅改變了裂紋尖端區(qū)域應(yīng)力場(chǎng),而且改變了梁結(jié)構(gòu)動(dòng)力特性.這些微觀裂紋若不及時(shí)發(fā)現(xiàn)、處理,其失穩(wěn)擴(kuò)展而引起的斷裂,將導(dǎo)致整個(gè)梁結(jié)構(gòu)產(chǎn)生毀滅性的破壞.故本文將針對(duì)移動(dòng)質(zhì)量-梁系統(tǒng)分析不同程度裂紋梁的振動(dòng)特性,對(duì)更加準(zhǔn)確地獲得梁結(jié)構(gòu)的固有頻率、保證梁結(jié)構(gòu)安全可靠具有重要意義.

對(duì)于移動(dòng)質(zhì)量-梁模型的研究,在過(guò)去的十幾年里研究人員在系統(tǒng)建模、動(dòng)力學(xué)分析方面已經(jīng)有許多研究成果[1-10],其中Wayou[7]研究了歐拉-伯努利梁在移動(dòng)力下的非線性動(dòng)態(tài)分析,彭獻(xiàn)[8-9]分析了移動(dòng)質(zhì)量作用下梁固有頻率的變化,Ismail Gerdemeli[10]利用移動(dòng)有限元法研究了隨著移動(dòng)質(zhì)量大小的變化對(duì)梁固有頻率的影響.

如何有效地識(shí)別和檢測(cè)裂紋具有十分重要的工程意義[11].對(duì)此,國(guó)內(nèi)外學(xué)者應(yīng)用不同的方法對(duì)有關(guān)裂紋梁的動(dòng)力特性問(wèn)題進(jìn)行了深入研究,從解析解、半解析解到數(shù)值方面的計(jì)算都提出了一系列的解決方案[13-20].而基于振動(dòng)特性的裂紋檢測(cè)法作為一種無(wú)損檢測(cè)的方法近年來(lái)越來(lái)越受到重視[12].

本文基于移動(dòng)有限元法,根據(jù)移動(dòng)質(zhì)量-梁系統(tǒng)的動(dòng)力學(xué)方程,建立了Euler-Bernoulli梁與移動(dòng)質(zhì)量耦合單元的時(shí)變質(zhì)量、剛度、阻尼矩陣.并利用卡式定理導(dǎo)出含裂紋梁?jiǎn)卧娜岫染仃?基于力平衡條件導(dǎo)出裂紋梁?jiǎn)卧膭偠染仃?在此基礎(chǔ)上,以堆場(chǎng)內(nèi)集裝箱起重機(jī)為例,計(jì)算了移動(dòng)質(zhì)量作用下含裂紋簡(jiǎn)支梁的動(dòng)態(tài)特性及動(dòng)力響應(yīng),并對(duì)結(jié)果進(jìn)行了對(duì)比分析.對(duì)集裝箱起重機(jī)的動(dòng)態(tài)設(shè)計(jì)、動(dòng)力性能評(píng)估以及振動(dòng)的控制具有一定的指導(dǎo)意義.

1 裂紋梁振動(dòng)的移動(dòng)有限元法

如圖1所示,將梁分為若干個(gè)單元,移動(dòng)質(zhì)量mp以恒定的速度vm沿梁運(yùn)動(dòng),其中第S個(gè)單元含有移動(dòng)質(zhì)量,第K個(gè)單元含有裂紋,考慮移動(dòng)慣性載荷作用下的系統(tǒng)運(yùn)動(dòng)方程為：

(1)

圖1 移動(dòng)質(zhì)量下含有裂紋的簡(jiǎn)支梁Fig. 1 Simple support crack beam due to a moving mass

1.1 移動(dòng)有限元原理

圖2為一載有移動(dòng)質(zhì)量mp的梁?jiǎn)卧?單元有兩個(gè)節(jié)點(diǎn),每個(gè)節(jié)點(diǎn)有兩個(gè)自由度.移動(dòng)質(zhì)量在梁上的位置隨著時(shí)間的變化而變化.

圖2 含有移動(dòng)質(zhì)量的梁?jiǎn)卧狥ig. 2 Beam element with moving mass

當(dāng)梁在振動(dòng)時(shí),梁在移動(dòng)質(zhì)量下的所受的橫向力為[3]：

(2)

其中：

(3)

fy(x,t)是移動(dòng)質(zhì)量在t時(shí)刻,x位置處作用的力.δ(x-xp)為狄拉克函數(shù),g為重力加速度.x0和v0是初始位置和初始速度.

考慮移動(dòng)質(zhì)量的慣性影響,梁的橫向加速度方程可以寫(xiě)為[19]：

(4)

當(dāng)質(zhì)量塊勻速運(yùn)動(dòng)時(shí),(4)式可以寫(xiě)為：

(5)

(5)式的另外一種形式為：

(6)

其中：“ .”為對(duì)時(shí)間的導(dǎo)數(shù),“ ′”為對(duì)位移的導(dǎo)數(shù),wy=wy(x,t)是梁的縱向撓度關(guān)于時(shí)間和位移的函數(shù).將(6)帶入(2)中,得：

(7)

在質(zhì)量塊勻速移動(dòng)下,單元的等效節(jié)點(diǎn)力為：

(8)

Ni(i=1～4)為梁?jiǎn)卧男魏瘮?shù)：

N1=1-3ξ(t)2+2ξ(t)3

N2=[ξ(t)-2ξ(t)2+ξ(t)3]l

N3=3ξ(t)2-2ξ(t)3

N4=[-ξ(t)2+ξ(t)3]l

(9)

形函數(shù)和位移、時(shí)間之間的關(guān)系為：

wy(x,t)=N1us1+N2us2+N3us3+N4us4

(10)

其中,usi(i=1～4)為單元的節(jié)點(diǎn)位移.

將(10)帶入(8),并將所得表達(dá)式寫(xiě)為矩陣形式：

(11)

其中,

(12e)

(12f)

(12g)

[m]、[c]、[k]為移動(dòng)有限元的時(shí)變質(zhì)量、阻尼和剛度矩陣.移動(dòng)質(zhì)量的位置xp(t)與質(zhì)量和速度有關(guān).

1.2 裂紋梁?jiǎn)卧?/p>

利用局部柔度理論建立裂紋梁?jiǎn)卧挠邢拊Ｐ?對(duì)于一個(gè)含有裂紋的歐拉-伯努利梁?jiǎn)卧?單元長(zhǎng)度為l,相對(duì)于梁?jiǎn)卧淖蠖肆鸭y位置為l1.應(yīng)用局部能量法來(lái)推導(dǎo)裂紋梁?jiǎn)卧膭偠染仃嘯20].

對(duì)于無(wú)裂紋梁?jiǎn)卧膽?yīng)變能為：

(13)

其中：E為楊氏模量,I為截面的慣性矩,l為梁?jiǎn)卧拈L(zhǎng)度.梁?jiǎn)卧膹澗胤匠炭梢詫?xiě)為：

Mb=Tl+M

(14)

其中：T為單元所受的剪力,M為單元的彎矩.

由于裂紋的存在,裂紋梁?jiǎn)卧鄬?duì)于一般的梁?jiǎn)卧嬖诟郊拥膽?yīng)變能,根據(jù)Tada[17]的理論：

(15)

其中：A為梁?jiǎn)卧臋M截面積,μ為泊松比,KⅠ、KⅡ、KⅢ為應(yīng)力強(qiáng)度因子.

圖3 裂紋梁?jiǎn)卧狥ig. 3 Cracked beam element

根據(jù)無(wú)裂紋梁理論,無(wú)裂紋梁?jiǎn)卧膹椥韵禂?shù)可以寫(xiě)為：

(16)

而考慮裂紋梁的附加應(yīng)變能,其彈性系數(shù)為：

(17)

對(duì)于一個(gè)裂紋梁?jiǎn)卧?單元的平衡方程可以寫(xiě)為：

(18)

根據(jù)公式(13)(14)和(15)～(17)彈性矩陣[C]可以寫(xiě)為：

(19)

根據(jù)公式(19),裂紋梁的剛度矩陣可以寫(xiě)為：

(20)

2 系統(tǒng)動(dòng)力方程求解

由圖1所示的多自由度系統(tǒng)的振動(dòng)方程為：

(21)

(22)

(23)

除了含有移動(dòng)質(zhì)量的S單元和含有裂紋的K單元：

(24)

(25)

(26)

(27)

隨時(shí)間變化的xm(t)和S可由(28)(29)得來(lái)：

xm(t)=xp(t)-(s-1)l

(28)

(29)

(30)

對(duì)于公式(21)所給出的動(dòng)力方程,可以通過(guò)Newmark數(shù)值積分來(lái)得出結(jié)果.

3 數(shù)值結(jié)果

本文給出的所有結(jié)果,是在重力加速度為9.81m/s2、阻尼率為ξ1=ξ2=0.05對(duì)應(yīng)于以前兩階固有頻率ω1,ω2；取m為移動(dòng)質(zhì)量的大小,M為梁的質(zhì)量.取梁的參數(shù)為：高h(yuǎn)=0.5m,寬b=1m,長(zhǎng)L=30m.簡(jiǎn)支梁的材料屬性：密度7890kg/m3、楊氏模量2.01×1011Pa、泊松比0.3.

3.1 移動(dòng)質(zhì)量下裂紋梁的固有頻率

為了揭示耦合系統(tǒng)各階頻率隨移動(dòng)質(zhì)量與梁的質(zhì)量比、位置比和裂紋相對(duì)位置的變化規(guī)律,按特征值計(jì)算法的結(jié)果,繪出該耦合系統(tǒng)的前三階固有頻率隨裂紋位置比β和系統(tǒng)質(zhì)量位置比λ變化的三維圖形分別如圖4、圖5、圖6所示.

圖4 耦合系統(tǒng)的第一階固有頻率隨λ和β變化的三維圖Fig. 4 3-dimension graph of 1st natural frequency of the coupled system with the change of λ and β

圖5 耦合系統(tǒng)的第二階固有頻率隨λ和β變化的三維圖Fig. 5 3-dimension graph of 2nd natural frequency of the coupled system with the change of λ and β

圖6 耦合系統(tǒng)的第三階固有頻率隨λ和β變化的三維圖Fig. 6 3-dimension graph of 3rd natural frequency of the coupled system with the change of λ and β

假定當(dāng)相對(duì)裂紋深度為0.5，質(zhì)量比為:m/M=0.2時(shí).由以下三圖可以清晰地看到,移動(dòng)質(zhì)量與簡(jiǎn)支梁耦合系統(tǒng)的各階固有頻率隨裂紋位置比和質(zhì)量位置比呈簡(jiǎn)諧函數(shù)規(guī)律變化,一階為半波,二階為全波,三階為1.5倍全波,且階數(shù)愈高波數(shù)越高.且根據(jù)圖像數(shù)據(jù)可以計(jì)算出,第一階頻率的變化范圍為0～25.49%；第二階頻率的變化范圍為0～17.31%；第三階頻率的變化范圍為0～14.71%.且當(dāng)裂紋位置和移動(dòng)質(zhì)量位置重合時(shí)的固有頻率相對(duì)于其他位置時(shí)會(huì)有較為顯著的下降.

綜上所述,移動(dòng)質(zhì)量與裂紋梁耦合系統(tǒng)的各階固有頻率并非常數(shù),而與系統(tǒng)的質(zhì)量比、位置比和裂紋相對(duì)深度、相對(duì)位置有關(guān)；若用梁的固有頻率代替耦合系統(tǒng)的固有頻率有時(shí)會(huì)產(chǎn)生較大的誤差.

3.2 移動(dòng)質(zhì)量下裂紋梁的動(dòng)態(tài)響應(yīng)

利用上述理論編寫(xiě)程序,將簡(jiǎn)支梁分為50個(gè)單元.取裂紋相對(duì)深度為0.5,相對(duì)位置為0.5,移動(dòng)質(zhì)量與梁的質(zhì)量比為：m/M=0.2.

圖7為移動(dòng)質(zhì)量在不同速度時(shí),梁跨中處的動(dòng)態(tài)響應(yīng).當(dāng)移動(dòng)質(zhì)量在裂紋梁上運(yùn)行時(shí),梁的撓度相比于無(wú)損梁的撓度加大.并且在同一移動(dòng)速度下,裂紋梁跨中處達(dá)到最大撓度的時(shí)刻要比無(wú)損梁達(dá)到最大撓度的時(shí)刻推后.并且隨著速度的增加,裂紋梁跨中處達(dá)到最大撓度的時(shí)刻有著向后推移的趨勢(shì).當(dāng)v=4m/s時(shí),移動(dòng)力與移動(dòng)質(zhì)量的計(jì)算誤差為6.4%；當(dāng)v=10m/s時(shí),移動(dòng)力與移動(dòng)質(zhì)量的計(jì)算誤差為18.5%.表明隨著速度的增加,移動(dòng)質(zhì)量慣性力對(duì)梁動(dòng)態(tài)響應(yīng)的影響增加.

圖7 動(dòng)質(zhì)量在不同速度作用下裂紋梁的動(dòng)態(tài)響應(yīng)(裂紋相對(duì)深度α=0.5；相對(duì)位置β=0.5)(--移動(dòng)力無(wú)裂紋, —○—移動(dòng)力有裂紋—移動(dòng)質(zhì)量無(wú)裂紋,—●—移動(dòng)質(zhì)量有裂紋)Fig. 7 Dynamic response due to moving mass in different velocity

圖8表明裂紋對(duì)梁在不同時(shí)刻振動(dòng)形狀的影響(v=10m/s、m/M=0.4).當(dāng)移動(dòng)質(zhì)量接近梁的右端時(shí),裂紋對(duì)梁振動(dòng)的影響越來(lái)越小.隨著質(zhì)量比和速度的增加,裂紋對(duì)梁撓屈變形影響最大的時(shí)刻并不是質(zhì)量移動(dòng)到裂紋位置時(shí)(裂紋相對(duì)位置為0.5),而是在t/T=0.6左右.

圖8 移動(dòng)質(zhì)量在速度為10m/s時(shí),裂紋梁和無(wú)裂紋梁在七個(gè)時(shí)刻的形狀Fig. 8 Beam shape for v=10m/s, α=0.5,mid-span

圖9所示為裂紋程度的加深對(duì)梁跨中處動(dòng)態(tài)響應(yīng)的影響(v=10m/s、m/M=0.4).隨著裂紋程度的增加跨中處動(dòng)態(tài)響應(yīng)的最大位移增大,當(dāng)裂紋相對(duì)深度α=0.2為時(shí),跨中處最大撓度值較無(wú)損梁下降了0.16%；當(dāng)α=0.3時(shí),跨中處最大撓度值較無(wú)損梁下降了0.64%；當(dāng)α=0.4時(shí),跨中處最大撓度值較無(wú)損梁下降了3.23%；當(dāng)α=0.5時(shí),跨中處最大撓度值較無(wú)損梁下降了13.09%.表明隨著裂紋相對(duì)深度的增加,裂紋對(duì)梁跨中處的動(dòng)態(tài)響應(yīng)越敏感.

圖9 裂紋深度的變化對(duì)動(dòng)態(tài)響應(yīng)的影響Fig. 9 Effect of crack size on the time response at mid-span

4 結(jié)論

本文利用移動(dòng)有限元和裂紋梁?jiǎn)卧碚搶?duì)移動(dòng)質(zhì)量下簡(jiǎn)支裂紋梁進(jìn)行動(dòng)力學(xué)計(jì)算分析,得出移動(dòng)質(zhì)量作用下裂紋的存在相比于無(wú)損梁會(huì)產(chǎn)生更大的動(dòng)態(tài)撓屈變形.運(yùn)用移動(dòng)有限元法,考慮移動(dòng)質(zhì)量的慣性力對(duì)梁結(jié)構(gòu)動(dòng)態(tài)性能的影響,隨著動(dòng)載速度的增加,相比于移動(dòng)力,慣性力對(duì)梁動(dòng)態(tài)特性的影響更顯著,且慣性力對(duì)裂紋的影響使得梁跨中處達(dá)到最大撓度的時(shí)刻向后偏置.當(dāng)裂紋相對(duì)深度達(dá)到α=0.5時(shí),梁的動(dòng)態(tài)撓度會(huì)有顯著的增大,若不及時(shí)發(fā)現(xiàn)處理將會(huì)導(dǎo)致結(jié)構(gòu)的破壞.

1 李海斌,畢世華,方遠(yuǎn)喬. 運(yùn)動(dòng)剛體激勵(lì)下彈性梁的振動(dòng)響應(yīng)分析. 工程力學(xué), 1998,15(4):74～80 (Li H B, Bi S H. Fang Y Q. Vibraction analysis of beams activated by moving rigid body.EngineeringMechanics, 1998,15(4):74～80 (in Chinese))

2 M.A.M Ahmoud, M.A.A Bou, Z Aid. Dynamic response of a beam with a cracksubject to a moving mass.JournalofSoundandVibration, 2002,256(4):591～603

3 Cifuentes A O. Dynamic response of a beam excited by a moving mass.FiniteElementsinAnalysisandDesign, 1989,5:237～246

4 Wu J J, Whittaker A R, Cartmell M P. The use of finite element techniques for calculating the dynamic response of structures to moving loads.ComputersandStructures, 2000:78:789～799

5 Foda M A, Abduljabbar Z. A dynamic Green's function formulation for the response of a beam structure to a moving mass.JournalofSoundandVibration, 1998,210(3):295～3066 Zhu X Q, Law S S. Precise time-step integration for the dynamic response of acontinuous beam under moving loads.JournalofSoundandVibration, 2001,240(5):962～970

7 Wayou A N Y, Tchoukuegno R, Woafo P. Non-linear dynamics of an elastic beam under moving loads.JournalofSoundandVibration, 2004,273:1101～1108

8 彭獻(xiàn),游福賀. 移動(dòng)質(zhì)量與梁耦合系統(tǒng)固有頻率的計(jì)算與分析. 動(dòng)力學(xué)與控制學(xué)報(bào), 2009,7(3):270～274 (Peng X, You F H. Calculation and analysis on natural frequency of a moving mass and beam′s coupled system.JournalofDynamicsandControl, 2009,7(3):270～274(in Chinese))

9 彭獻(xiàn),殷新鋒. 車-橋系統(tǒng)的振動(dòng)分析及控制. 動(dòng)力學(xué)與控制學(xué)報(bào), 2006,4(3):253～258 (Peng X, Yin X F. Vibration analysis and control on vehicle-bridge system.JournalofDynamicsandControl, 2006,4(3):253～258)

10Ismail Gerdemeli, Ismail Esen, Derya?zer. Dynamic response of an overhead crane beam due to a moving mass using moving finite element approximation.KeyEngineeringMaterials, 2011,450:99～102

11范天佑. 斷裂理論基礎(chǔ). 北京：科學(xué)出版社, 2003 (Fan T Y. Theoretical basis of fracture. Beijing: Press of Science, 2003 (in Chinese))

12吳寧祥,謝里陽(yáng),吳克勤,由美雁. 裂紋梁動(dòng)力學(xué)仿真的有限元模型. 東北大學(xué)學(xué)報(bào), 2008,29(5):714～717 (Wu N X, Xie L Y, Wu K Q, You M Y. An improved finite element model for dynamic simulation ofcracked beam.JournalofNortheasternUniversity, 2008,29(5)714～717(in Chinese))

13劉文光,陳國(guó)平. 呼吸式裂紋梁的振動(dòng)疲勞裂紋擴(kuò)展耦合分析. 中國(guó)機(jī)械工程, 2010,21(23):2798～2802 (Liu W G, Chen G P.Coupling analysis of vibration fatigue crack growth for breathing cracked beam.MechanicalEngineeringofChina, 2010,21(23):2798～2802(in Chinese))14李兆軍,龍慧,劉洋,邱旻. 基于有限元位移模式的含裂紋梁結(jié)構(gòu)動(dòng)力學(xué)模型. 中國(guó)機(jī)械工程, 2014,25(12):1563～1566 (Li Z J, Long H,Liu Y, Qiu H. Dynamics equation of cracked beam based on finite element displacement mode.MechanicalEngineeringofChina,2014,25(12):1563～1566(in Chinese))

15吳寧祥,謝里陽(yáng),吳克勤等. 裂紋梁的相似特性及模型誤差分析. 機(jī)械強(qiáng)度, 2006,28(4):570～573 (Wu N X, Xie L Y, Wu K Q, et al. Similarity of cracked beam and model error analysis.MechanicalStrength, 2006,28(4):570～573 (in Chinese))

16錢(qián)管良,顧松年,姜節(jié)勝. 含裂紋梁的動(dòng)力響應(yīng). 振動(dòng)工程學(xué)報(bào), 1989,2(3):78～85 (Qian G K, Gu S N, Jiang J S. The dynamic response of beam with crack.JournalofVibrationEngineering, 1989,2(3):78～85 (in Chinese))

17Tada H, Paris P, Irwin G. The stress analysis of cracks handbook. Hellertown, PA: Del Research Corp, 1973

18Clough R W, Penzien J. Dynamics of structures 3rd ed. Computers and Structures, Inc., Berkeley, 2003

19Fryba L. Vibration solids and structures under moving loads, Thomas TelfordHouse, London, 1999

20Nahvi H . Crack detection in beams using experimental modal data and finite element model. Mechanical Science, 2000

? Corresponding author E-mail:dick77891@qq.com

23 January 2016,revised 10 October 2016.

VIBRATION ANALYSIS OF A CRACKED BEAM BASES ON MOVING FINITE ELEMENT APPROACH*

Liang Gang1?Cao Pei1Tang Jiong2Shen Qing2,Li Xinshuang2

(1.CollegeofLogisticsEngineering,ShanghaiMaritimeUniversity,Shanghai201399,China)(2.ShanghaiZhenhuaHeavyIndustryCoLtd,Shanghai200125,China)(3.DepartmentofMechanicalEngineering,UniversityofConnecticut,American)

The vibration of a cracked beam due to a moving mass is investigated by using the moving finite element method and the local flexibility method. The effect of the relative position and relative depth of the crack on the natural frequency and the dynamic response of the cracked beam under different velocity are discussed through the numerical calculation. The analysis shows that the crack and the moving mass increase the dynamic displacement of beam, but decrease the coupling frequency of the system in a certain degree. Meanwhile, with the increase of the speed of moving mass and the relative depth of crack, the dynamic response of the beam is apt to increase.

moving finite element, simple support beam, open crack, natural frequency, dynamic response

10.6052/1672-6553-2016-051

2016-01-23收到第1稿,2016-10-10收到修改稿.

? 通訊作者 E-mail:dick77891@qq.com