任慧龍，于鵬垚，李 輝，馮國(guó)慶

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001）

船體三維變形響應(yīng)的數(shù)值預(yù)報(bào)

任慧龍，于鵬垚，李 輝，馮國(guó)慶

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001）

針對(duì)航行中艦船變形問(wèn)題，采用三維水彈性理論實(shí)現(xiàn)了波浪中艦船的動(dòng)態(tài)變形預(yù)報(bào)。通過(guò)擬合不同模態(tài)下船體有限元模型的節(jié)點(diǎn)位移，實(shí)現(xiàn)了將三維模態(tài)信息應(yīng)用于廣義流固耦合界面條件。建立了艦船在規(guī)則波中的廣義水彈性運(yùn)動(dòng)方程，求解出各階振動(dòng)模態(tài)的主坐標(biāo)。將各階彈性模態(tài)主坐標(biāo)與預(yù)報(bào)位置的角變形矩陣相乘，得到角變形響應(yīng)。針對(duì)某實(shí)船，分別采用三維水彈性方法和三維剛體方法進(jìn)行計(jì)算，通過(guò)對(duì)船體結(jié)構(gòu)角變形響應(yīng)的對(duì)比，驗(yàn)證了三維水彈性方法的正確性，同時(shí)也證明了該方法在計(jì)算效率上的優(yōu)勢(shì)。

艦船變形；流固耦合；三維水彈性方法；模態(tài)分析；變形響應(yīng)；數(shù)值預(yù)報(bào)

隨著船舶的大型化，船體在波浪中航行時(shí)，由周圍波動(dòng)壓力作用引起的彈性變形愈發(fā)明顯。對(duì)于某些特殊用途的船舶（如科考船、軍艦等），較大的動(dòng)態(tài)變形會(huì)影響設(shè)備的正常使用，進(jìn)而影響其正常作業(yè)。因此，國(guó)內(nèi)外相關(guān)學(xué)者提出一些方法對(duì)船舶營(yíng)運(yùn)時(shí)的波浪載荷誘導(dǎo)的動(dòng)態(tài)變形響應(yīng)進(jìn)行預(yù)報(bào)。陳曙梅等［1?4］分2步進(jìn)行載荷預(yù)報(bào)和結(jié)構(gòu)響應(yīng)計(jì)算，這些學(xué)者先將艦船看作三維剛體結(jié)構(gòu)來(lái)計(jì)算不同工況下波浪載荷，再將載荷加載到三維有限元結(jié)構(gòu)上來(lái)獲取船體變形，這種方法在后續(xù)論述中簡(jiǎn)稱為剛體方法。陳敏康等［5?6］等采用二維水彈性理論對(duì)船體的變形進(jìn)行預(yù)報(bào)。

剛體方法是在準(zhǔn)靜態(tài)假定下進(jìn)行計(jì)算，無(wú)法反映真實(shí)船體結(jié)構(gòu)動(dòng)力響應(yīng)特性，而且在計(jì)算船體角變形幅頻響應(yīng)時(shí)，需要針對(duì)較多頻率進(jìn)行壓力加載和靜力分析，工作量較大。二維水彈性理論是將船體結(jié)構(gòu)簡(jiǎn)化為一根變截面的船體梁，不能夠考慮船寬方向上船體結(jié)構(gòu)的變形。本文利用三維有限元模型的模態(tài)分析結(jié)果，采用的三維水彈性預(yù)報(bào)方法實(shí)現(xiàn)了船體三維變形響應(yīng)的數(shù)值預(yù)報(bào)。

1 船體角變形矩陣的建立

對(duì)角變形響應(yīng)而言，關(guān)心的是一個(gè)位置相對(duì)另一位置的角變形，如船體上某位置相對(duì)于羅經(jīng)平臺(tái)的變形。本文引入角變形矩陣X3n×m來(lái)表達(dá)不同模態(tài)下不同位置間的相對(duì)角變形振型。為了便于進(jìn)行描述，引入如圖1所示坐標(biāo)系，其中原點(diǎn)位于重心在靜水面上的投影，x軸有船舵指向船艏，y由右舷指向左舷，z軸數(shù)值向上。

圖1 坐標(biāo)系Fig.1 Coordinate system

本文采用三維有限元方法對(duì)船體結(jié)構(gòu)進(jìn)行干模態(tài)分析，并以艉垂線與上甲板交點(diǎn)的線位移進(jìn)行歸一化，得到各階彈性模態(tài)對(duì)應(yīng)的振型，廣義質(zhì)量矩陣a，廣義剛度矩陣c。

2 流固耦合界面條件的數(shù)值處理

基于勢(shì)流理論求解艦船周圍的流場(chǎng)時(shí)，需要求解流場(chǎng)的繞射勢(shì)φd和輻射勢(shì)φk。通?？刹捎迷磪R分布法［7?9］，結(jié)合自由面條件、遠(yuǎn)方條件和廣義的流固耦合的界面條件［10］進(jìn)行求取。

廣義的流固耦合界面條件：

式中：φ0為已知的單位波幅入射勢(shì)；ω為波浪自然頻率；n為面元單位法向量；uk為第k模態(tài)下的面元中心點(diǎn)位移；W表示船體周圍流場(chǎng)的定常速度，忽略定常速度勢(shì)時(shí)，W＝-Ui，i為x坐標(biāo)軸的單位方向向量，U為航速。

從式（2）可以看出，船體結(jié)構(gòu)的振型信息是通過(guò)面元中心點(diǎn)的位移和該位移對(duì)x、y、z方向的偏導(dǎo)數(shù)體現(xiàn)在彈性結(jié)構(gòu)的水動(dòng)力求解中。因而，問(wèn)題的關(guān)鍵就是如何通過(guò)三維模型在不同模態(tài)下的節(jié)點(diǎn)的位移信息來(lái)獲得水動(dòng)力求解中的面元中心點(diǎn)的位移和位移對(duì)x、y、z方向的偏導(dǎo)數(shù)。本文采用擬合位移形函數(shù)的方法來(lái)解決這一問(wèn)題。

船體結(jié)構(gòu)三維有限元模型的節(jié)點(diǎn)線位移矢量u＝（ux，uy，uz）包含3個(gè)方向的位移分量，本文以u(píng)x方向位移為例，推導(dǎo)出如何利用三維有限元模型外殼的振型信息確定對(duì)應(yīng)模態(tài)下面元中心點(diǎn)的位移及其導(dǎo)數(shù)，其他兩個(gè)方向位移的推導(dǎo)方法與其相同，不再贅述。

為便于后續(xù)描述，作如下定義：將面元中心點(diǎn)到面元節(jié)點(diǎn)的最遠(yuǎn)距離定義為面元控制域半徑，以該半徑作球的區(qū)域定義為面元控制域。則面元中心點(diǎn)的位移和位移的導(dǎo)數(shù)由面元控制域內(nèi)節(jié)點(diǎn)的振型信息來(lái)確定。

對(duì)于空間運(yùn)動(dòng)的剛體，當(dāng)給定剛體質(zhì)心處的6自由度運(yùn)動(dòng)，便可以得到任意位置的位移：

式中：un為待求任意一點(diǎn)的線位移；uc、ωc分別為質(zhì)心處的線位移和角位移；r為待求位置與質(zhì)心之間的向徑；xn、yn、zn為待求位置的坐標(biāo)；xc、yc、zc為質(zhì)心坐標(biāo)；i、j、k分別為x、y、z坐標(biāo)軸的單位方向向量。

可見(jiàn)，任意節(jié)點(diǎn)線位移可由節(jié)點(diǎn)坐標(biāo)與剛體質(zhì)心處運(yùn)動(dòng)相乘的形式來(lái)表示。借鑒這種思路，當(dāng)面元不與坐標(biāo)軸垂直時(shí)，假定船體外殼上，面元對(duì)應(yīng)區(qū)域內(nèi)不同位置的x方向位移形函數(shù)具有如下形式：

式中：a0、a1、a2、a3為待定系數(shù)。

當(dāng)面元與坐標(biāo)軸垂直時(shí)，如與z軸垂直，控制域內(nèi)的各個(gè)節(jié)點(diǎn)可能為同一個(gè)值，若同時(shí)為零，則a3為任意值，可導(dǎo)致后續(xù)求解方程組時(shí)系數(shù)矩陣奇異。故當(dāng)面元與坐標(biāo)軸垂直時(shí)，假定型函數(shù)的形式如下：

可以看出想要確定位移型函數(shù)，面元控制域內(nèi)節(jié)點(diǎn)不得少于4，當(dāng)少于4時(shí)可以適當(dāng)增加控制域半徑來(lái)增加節(jié)點(diǎn)。

若控制域內(nèi)包含n個(gè)節(jié)點(diǎn)，可根據(jù)振型信息和位移型函數(shù)確定n個(gè)方程。方程的個(gè)數(shù)往往大于待定系數(shù)，即n＞4，此時(shí)可以基于最小二乘理論確定方程的系數(shù)，進(jìn)而確定型函數(shù)，代入面元中心點(diǎn)坐標(biāo)得到面元中心點(diǎn)位移。型函數(shù)分別對(duì)x、y、z方向求導(dǎo)，代入面元中心點(diǎn)坐標(biāo)便得到面元中心點(diǎn)位移的導(dǎo)數(shù)。

對(duì)于與坐標(biāo)軸垂直的面元，這里定性的確定該面元中心點(diǎn)位移對(duì)所垂直坐標(biāo)軸方向的導(dǎo)數(shù)。如當(dāng)面元垂直x軸時(shí)，該面元應(yīng)位于結(jié)構(gòu)的端部，通常對(duì)應(yīng)于船體艏艉兩端，該處對(duì)應(yīng)的剖面內(nèi)力為零，則各方向位移對(duì)x的導(dǎo)數(shù)應(yīng)為零。同理對(duì)于垂直于y軸或z軸的面元，各方向位移對(duì)y軸或z軸的導(dǎo)數(shù)應(yīng)為零。

3 船體角變形響應(yīng)的求解

根據(jù)水彈性力學(xué)建立起適用于規(guī)則波激勵(lì)下的廣義線性水彈性運(yùn)動(dòng)方程［11］為

式中：a為船體廣義質(zhì)量矩陣，b為船體廣義剛度矩陣，a和c可由前面的三維有限元結(jié)構(gòu)模態(tài)分析獲得；b為船體廣義阻尼矩陣，可通過(guò)試驗(yàn)或經(jīng)驗(yàn)公式獲得；pa為主坐標(biāo)幅值列向量；ωe為遭遇頻率；A、B、C為廣義水動(dòng)力系數(shù)，具體含義分別為廣義流體附加質(zhì)量、廣義流體附加阻尼和廣義流體恢復(fù)力系數(shù)矩陣；F為廣義波浪激勵(lì)力。

對(duì)于廣義水動(dòng)力系數(shù)A、B、C和廣義波浪激勵(lì)力F，具體可依據(jù)下式求得：

式中：ρ為流體的密度；ρb為船體結(jié)構(gòu)材料的密度；ds為船體濕表面的面積微元；dV為船體結(jié)構(gòu)的體積微元；n為面元中心點(diǎn)法向量；ur為第r階模態(tài)下的面元中心點(diǎn)位移；φk為第k階模態(tài)對(duì)應(yīng)的輻射勢(shì)；wk為第k模態(tài)下面元中心點(diǎn)的垂向位移；gs表示重力加速度矢量，gs＝（0，0，-g）；θk表示第k階船體振型產(chǎn)生的結(jié)構(gòu)轉(zhuǎn)角變形，θk＝（?×uk）／2；ζa為入射波波幅；φ0為已知的單位波幅入射勢(shì)；φd為單位波幅繞射勢(shì)；U為航速。

綜上，在確定方程（7）的各項(xiàng)系數(shù)后，可求得在波浪激勵(lì)下不同模態(tài)的主坐標(biāo)響應(yīng)。主坐標(biāo)響應(yīng)與角變形矩陣X3n×m相乘便可得到船體的任意位置間的角變形響應(yīng)。

4 應(yīng)用實(shí)例

利用提出的基于三維水彈性理論預(yù)報(bào)變形的方法，對(duì)某艦船在波浪作用下的角變形響應(yīng)進(jìn)行預(yù)報(bào)。首先，將三維水彈性得到的運(yùn)動(dòng)響應(yīng)與傳統(tǒng)的三維剛體載荷方法進(jìn)行比較，然后再將剛體方法得到的波浪載荷加載到三維有限元模型上進(jìn)行準(zhǔn)靜態(tài)分析，并提出取船體上任意位置間的相對(duì)角變形響應(yīng)，最后將兩種方法得到的角變形響應(yīng)進(jìn)行比較。

考慮到對(duì)于船體上任意位置間的相對(duì)角變形響應(yīng)預(yù)報(bào)方法相同，本文僅以船體甲板上某一位置相對(duì)于船舯附近的羅經(jīng)平臺(tái)的變形為例進(jìn)行角變形響應(yīng)預(yù)報(bào)。為了便于設(shè)備的使用，預(yù)報(bào)點(diǎn)位置和羅經(jīng)平臺(tái)都在船體的強(qiáng)力構(gòu)件附近，并采取了局部加強(qiáng)。

4.1 艦船主尺度數(shù)據(jù)及關(guān)注位置

選取某船舶進(jìn)行船體三維變形響應(yīng)的數(shù)值預(yù)報(bào)的相關(guān)研究，船舶的主尺度數(shù)據(jù)和關(guān)注位置坐標(biāo)分別如表1、2所示。

表1 艦船主尺度數(shù)據(jù)Table1 Principal parameters of the ship

表2 關(guān)注位置坐標(biāo)Table2 Coordinates of concerned positions

4.2 干模態(tài)分析

本文采用有限元軟件MSC.nastran對(duì)船體整船結(jié)構(gòu)干模態(tài)進(jìn)行分析，考慮到關(guān)注的位置位于甲板上，在圖2中列出彈性模態(tài)對(duì)應(yīng)甲板位置的振型。表3列出各階模態(tài)固有頻率、廣義質(zhì)量和無(wú)因次阻尼系數(shù)。表4給出關(guān)注位置相對(duì)羅經(jīng)平臺(tái)的角變形矩陣。實(shí)際計(jì)算中，所選取的彈性模態(tài)的數(shù)目將在后續(xù)進(jìn)行進(jìn)一步的討論。對(duì)于船體自身的結(jié)構(gòu)阻尼，理論計(jì)算很難實(shí)現(xiàn)精確的確定，而且結(jié)構(gòu)阻尼主要影響船體與波浪共振時(shí)的響應(yīng)峰值，對(duì)低頻波浪激勵(lì)下的結(jié)構(gòu)響應(yīng)影響不大，因而本文計(jì)算中的阻尼系數(shù)主要參照相關(guān)文獻(xiàn)［12］的經(jīng)驗(yàn)公式進(jìn)行確定。本文僅以前4階彈性模態(tài)為例列出模態(tài)分析結(jié)果，對(duì)于更多的模態(tài)只需擴(kuò)展相應(yīng)的行或列。

圖2 彈性模態(tài)甲板振型Fig.2 Deck shape of elastic principal mode

表3 船體干模態(tài)固有頻率、廣義質(zhì)量和無(wú)因次阻尼系數(shù)Table3 Eigenfrequencies，generalized mass of flexible modes and nondimensional damping coefficients

表4 角變形矩陣Table4 Angular deformation matrixmin

由船體干模態(tài)分析結(jié)果可以看出，一階彈性模態(tài)振動(dòng)形式以兩節(jié)點(diǎn)的垂向振動(dòng)為主。若采用二維的船體梁模型進(jìn)行模態(tài)分析，對(duì)于對(duì)稱船體，該模態(tài)下任意位置處繞x和z軸的角變形均為0，而由三維有限元方法得到相對(duì)角變形矩陣可以看出，在一階彈性模態(tài)下，繞x和z軸的對(duì)應(yīng)元素并不0，說(shuō)明由于船體結(jié)構(gòu)的局部特點(diǎn)影響，在以垂向振動(dòng)為主模態(tài)下，對(duì)于繞x和z軸的角變形也有貢獻(xiàn)。

4.3 彈性模態(tài)數(shù)目的選取

對(duì)于實(shí)際的船體結(jié)構(gòu)，從理論上說(shuō)，它具有無(wú)限種振動(dòng)形式，然而在采用模態(tài)疊加方法計(jì)算某一種外力激勵(lì)下的結(jié)構(gòu)響應(yīng)時(shí)，不可能考慮到結(jié)構(gòu)所有模態(tài)，只能夠盡可能的選取對(duì)結(jié)構(gòu)響應(yīng)貢獻(xiàn)較大的模態(tài)用于計(jì)算。對(duì)于本文計(jì)算的船體結(jié)構(gòu)，結(jié)構(gòu)響應(yīng)主要是由于波浪的激勵(lì)引起的，而實(shí)際船體的固有頻率通常大于波浪頻率的分布范圍，因而固有頻率越低的模態(tài)越接近波浪頻率，對(duì)結(jié)構(gòu)響應(yīng)貢獻(xiàn)越大。本文為考察究竟多少個(gè)彈性模態(tài)能夠充分反映船體關(guān)注位置的角變形響應(yīng)，分別采用前4、8、12階彈性模態(tài)用于水彈性計(jì)算。計(jì)算時(shí)，艦船航速為20 kn，浪向角為30°。圖3為選取不同彈性模態(tài)數(shù)目進(jìn)行計(jì)算時(shí)的角變形幅頻響應(yīng)曲線。

通過(guò)計(jì)算可以看出，選取4個(gè)彈性模態(tài)數(shù)目的計(jì)算與選取8和12個(gè)彈性模態(tài)的計(jì)算結(jié)果差別較大，選取8個(gè)和12個(gè)彈性模態(tài)的計(jì)算結(jié)果幾乎相同，說(shuō)明選取8個(gè)彈性模態(tài)已經(jīng)能夠充分反映出該海況下船體關(guān)注位置的角變形響應(yīng)，其中第8階彈性模態(tài)對(duì)應(yīng)的干模態(tài)固有頻率為34 rad／s，大于波浪能夠激勵(lì)的結(jié)構(gòu)干模態(tài)頻率范圍。因而后續(xù)采用三維水彈性方法進(jìn)行計(jì)算時(shí)均采用8個(gè)彈性模態(tài)。

圖3 不同模態(tài)數(shù)目下角變形幅頻響應(yīng)曲線Fig.3 Angular deformation amplitude?frequency re?sponse curves of different mode numbers

4.4 變形響應(yīng)的對(duì)比

在航速20 kn，浪向角30°時(shí)，分別采用三維水彈性方法和加載三維剛體方法的波浪載荷的方法預(yù)報(bào)船體的角變形的幅頻響應(yīng)，如圖4。

圖4 2種方法的角變形幅頻響應(yīng)曲線Fig.4 Angular deformation amplitude?frequency re?sponse curves of the two methods

通過(guò)對(duì)比可以看出，2種計(jì)算結(jié)果的趨勢(shì)基本一致，而且繞y軸的角變形明顯大于繞x、z軸的角變形。對(duì)于繞x、z軸的角變形，剛體方法的計(jì)算結(jié)果大于水彈性方法，對(duì)于繞y軸的角變形，2種方法基本一致。分析如下：三維水彈性方法是一種考慮到流固耦合的動(dòng)力分析方法，而剛體方法是在船體是剛體，結(jié)構(gòu)響應(yīng)是準(zhǔn)靜態(tài)的假定下進(jìn)行的計(jì)算，因而2種方法的計(jì)算結(jié)果必然存在偏差?？傮w的看來(lái)，對(duì)于船體角變形比較大的方向（繞y軸角變形），2種方法預(yù)報(bào)結(jié)果幾乎相同，一定程度上驗(yàn)證了水彈性方法正確性。

5 結(jié)論

對(duì)于預(yù)報(bào)艦船在波浪中角變形這一流固耦合問(wèn)題，本文采用擬合位移形函數(shù)的處理方法，實(shí)現(xiàn)將三維有限元的模態(tài)分析結(jié)果應(yīng)于與三維水彈性理論來(lái)預(yù)報(bào)結(jié)構(gòu)的角變形響應(yīng)。通過(guò)計(jì)算表明：

1）本文所選用的前8階彈性模態(tài)已足夠可以反映出關(guān)注位置處的角變形響應(yīng)。

2）通過(guò)剛體方法和水彈性方法的變形響應(yīng)對(duì)比，驗(yàn)證了采用擬合位移形函數(shù)的方法處理流固耦合界面條件的正確性。

3）通過(guò)水彈性方法，避免了將水動(dòng)壓力加載到有限元模型上，可以快速的獲得艦船的變形幅頻響應(yīng)，方便了后續(xù)變形響應(yīng)的長(zhǎng)短期分析和可靠性分析。

［1］陳曙梅.航行狀態(tài)下艦船結(jié)構(gòu)總體變形分析［D］.哈爾濱：哈爾濱工程大學(xué)，2004：7?42.

CHEN Shumei.Deformation analysis of whole vessel in navi?gating state［D］.Harbin：Harbin Engineering University，2004：7?42.

［2］LIU Aili，MA Hongxu，Dai Hongde.Estimation method for ship deformation under wave Loads［C］／／International Con?ference on Mechanic Automation and Control Engineering.Wuhan，China，2010：93?96.

［3］MURAWSKI L.Shaft line alignment analysis taking ship construction flexibility and deformations into consideration［J］.Marine Structures，2005，18：62?84.

［4］YU Pengyao，F(xiàn)ENG Guoqing，REN Huilong，et al.Deform?ation analysis and reliability assessment of ship hull in irreg?ular waves［C］／／Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering.Nantes，F(xiàn)rances，2013.

［5］陳敏康，邱強(qiáng).隨機(jī)波浪中船體變形轉(zhuǎn)角的估算方法及應(yīng)用［J］.船舶力學(xué)，2006，10（2）：87?95.

CHEN Minkang，QIU Qiang.Estimate method and its appli?cations for angular deflection of girder in random wave［J］.Journal of Ship Mechanics，2006，10（2）：87?95.

［6］TAYLOR R E.Reflections on the structural dynamics of floating beams and ship hulls in waves［J］.Journal of Engi?neering for the Maritime Environment，2009，223（3）：257?273.

［7］劉應(yīng)中，繆國(guó)平.船舶在波浪中的運(yùn)動(dòng)理論［M］.上海：上海交通大學(xué)出版社，1987：56?80.

［8］戴遺山，段文洋.船舶在波浪中的勢(shì)流理論［M］.北京：國(guó)防工業(yè)出版社，2007：63?84.

［9］李輝.船舶波浪載荷的三維水彈性分析方法研究［D］.哈爾濱：哈爾濱工程大學(xué)，2009：46?66.

［10］夏錦祝，吳有生.流固耦合問(wèn)題的一個(gè)一般交界面條件［J］.艦船性能研究，1993（2）：15?22.

XIA Jinzhu，WU Yousheng.One common interface condi?tion of fluid solid interaction problem［J］.Ship Perform?ance Research，1993（2）：15?22.

［11］李輝，任慧龍，劉日明，等.FPSO運(yùn)動(dòng)與波浪載荷三維水彈性分析方法研究［J］.哈爾濱工程大學(xué)報(bào)，2008，29（8）：789?793.

LI Hui，REN Huilong，LIU Riming，et al.3?D hydroelas?tic method for analysis of motion and wave loads of an FP?SO ship［J］.Journal of Harbin Engineering University，2008，29（8）：789?793.

［12］董艷秋.船舶波浪外荷和水彈性［M］.天津：天津大學(xué)出版社，1991：104?141.

Numerical prediction of three dimensional deformation response of the ship hull

REN Huilong，YU Pengyao，LI Hui，F(xiàn)ENG Guoqing
（College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China）

Aiming at the deformation problem of the navigating ship，the three dimensional（3D）hydroelastic meth?od is adopted to predict the dynamic deformation of ships in waves.Through fitting the nodes displacement of differ?ent vibration modes，the modal information of the 3D finite model is applied to the generalized interface boundary condition of fluid?structure interactions.Then，generalized hydroelastic equations of ship motion are established and the principal coordinates of different order vibration modes are solved.Furthermore，the angular deformation re?sponse could be achieved through multiplying the angular deformation matrix and the principal coordinates of differ?ent order elastic vibration modes.Take a real ship as an example，the 3D hydroelastic method and the 3D rigid method are applied to calculate the deformation response of the ship，respectively.Through comparing the results from the two methods，the correctness and high efficiency of the 3D hydroelastic method are proved.

ship deformation；fluid and structural interactions；3D hydroelastic method；modal analysis；fluid?structure interactions；angular deformation matrix；deformation response；numerical prediction

10.3969／j.issn.1006?7043.201311053

U661.3

A

1006?7043（2015）01?0134?05

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006?7043.201311053.html

2013?11?15.網(wǎng)絡(luò)出版時(shí)間：2014?11?07.

國(guó)家自然科學(xué)基金資助項(xiàng)目（51079034）.

任慧龍（1965?），男，教授，博士生導(dǎo)師.

任慧龍，E?mail：renhuilong＠263.net.