船舶結(jié)構(gòu)B樣條小波無(wú)網(wǎng)格分析技術(shù)

陳建平，唐文勇，徐曼平

(1.上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海200240;2.廣州航海學(xué)院船舶工程學(xué)院，廣東廣州510725)

?

船舶結(jié)構(gòu)B樣條小波無(wú)網(wǎng)格分析技術(shù)

陳建平1，2，唐文勇1，徐曼平2

(1.上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海200240;2.廣州航海學(xué)院船舶工程學(xué)院，廣東廣州510725)

摘要:為了解決船舶平直結(jié)構(gòu)場(chǎng)量高梯度自適應(yīng)分析問(wèn)題，提出了基于B樣條小波的無(wú)網(wǎng)格局部Petrov－Galerkin法。首先運(yùn)用最小二乘法和加權(quán)余量法來(lái)求解結(jié)構(gòu)位移場(chǎng)量的逼近函數(shù)，并給出了問(wèn)題的控制方程和剛度方程。然后在局部無(wú)網(wǎng)格Petrov－Galerkin法的基礎(chǔ)上，利用m階B樣條函數(shù)作為小波基函數(shù)來(lái)構(gòu)造船舶結(jié)構(gòu)位移場(chǎng)的逼近函數(shù)，并采用兩尺度分解技術(shù)來(lái)分解應(yīng)力場(chǎng)的高梯度成分和低尺度成分，應(yīng)用高尺度成分來(lái)表示應(yīng)力高梯度成分。最后選取了兩種典型船舶結(jié)構(gòu)進(jìn)行變形和應(yīng)力分析，并通過(guò)與有限元法的計(jì)算結(jié)果進(jìn)行比較，驗(yàn)證了本文提出方法的有效性。

關(guān)鍵詞:船舶結(jié)構(gòu)分析;逼近函數(shù);移動(dòng)最小二乘法;B樣條小波;尺度函數(shù);無(wú)網(wǎng)格MLPG法

當(dāng)前船舶結(jié)構(gòu)分析最常用也是最為有效的計(jì)算工具之一就是有限元方法［1］。但有限元法因其自身的特點(diǎn)，在處理分析船舶結(jié)構(gòu)場(chǎng)量(位移場(chǎng)和應(yīng)力場(chǎng)等)變化劇烈的高梯度區(qū)域時(shí)，會(huì)出現(xiàn)計(jì)算精度降低甚至計(jì)算中斷現(xiàn)象，為了解決這個(gè)問(wèn)題，有限元法通常采用的方法就是對(duì)該區(qū)域的網(wǎng)格進(jìn)行加密(細(xì)分)或者采用高階單元，這樣就要求方法具有較強(qiáng)的自適應(yīng)分析能力，其結(jié)果就加大了有限元法前后處理的工作量，從而降低了計(jì)算效率，事實(shí)上這種方法也并不能從根本上消除問(wèn)題產(chǎn)生的根源。作為與有限元法相對(duì)應(yīng)的另一種數(shù)值分析方法——無(wú)網(wǎng)格法，在近20年中得到了很大的發(fā)展［2－6］。無(wú)網(wǎng)格法是建立在系列獨(dú)立的離散點(diǎn)的基礎(chǔ)上，通過(guò)構(gòu)造點(diǎn)的近似函數(shù)來(lái)求解問(wèn)題。與傳統(tǒng)的有限元法相比，無(wú)網(wǎng)格法無(wú)需網(wǎng)格背景，在計(jì)算過(guò)程中可以根據(jù)需要任意增減節(jié)點(diǎn)，而不需要處理節(jié)點(diǎn)之間的拓?fù)湫畔?，特別適合用來(lái)進(jìn)行自適應(yīng)分析計(jì)算。無(wú)網(wǎng)格法目前在航空材料、高速碰撞、動(dòng)態(tài)裂紋擴(kuò)展、加工成型、節(jié)理巖體分析等諸多領(lǐng)域都得到了較為廣泛的應(yīng)用［7－10］。無(wú)網(wǎng)格法的計(jì)算精度同有限元法一樣因其自身的特點(diǎn)，受節(jié)點(diǎn)密度、基函數(shù)的選用及其階次以及支撐域半徑大小等因素的影響［11－12］。文章用局部無(wú)網(wǎng)格Petrov－Galerkin法為基礎(chǔ)，利用m階B樣條函數(shù)作為小波基函數(shù)來(lái)構(gòu)造對(duì)船舶結(jié)構(gòu)應(yīng)力場(chǎng)的逼近函數(shù)，并采用兩尺度分解技術(shù)來(lái)分解應(yīng)力場(chǎng)的高梯度成分和低尺度成分。文章最后選取了兩種典型船舶結(jié)構(gòu)進(jìn)行變形和應(yīng)力分析，并與有限元法(ANSYS軟件)的計(jì)算結(jié)果進(jìn)行了比較來(lái)驗(yàn)證文章方法的有效性。

1　船舶結(jié)構(gòu)MLPG法基本表達(dá)

由于船體板的板厚遠(yuǎn)小于其另外兩個(gè)方向(長(zhǎng)度和寬度)的幾何尺寸，在分析其受力時(shí)，可以運(yùn)用Kirchhoff－Love板殼理論和Mindlin－Reissner板殼理論。文章基于板殼的Mindlin－Reissner基本理論，運(yùn)用局部無(wú)網(wǎng)格Petrov－Galerkin法(MLPG)來(lái)分析研究無(wú)網(wǎng)格法下的船舶板結(jié)構(gòu)應(yīng)力。

1.1節(jié)點(diǎn)控制方程

板結(jié)構(gòu)經(jīng)離散后二維彈性理論為基礎(chǔ)的節(jié)點(diǎn)系統(tǒng)方程［11］:

利用加權(quán)余量法［12］，可以得到節(jié)點(diǎn)I處系統(tǒng)方程的微分方程強(qiáng)形式為

對(duì)節(jié)點(diǎn)I依式(2)在其積分域上進(jìn)行積分，可以建立起它的系統(tǒng)方程。這樣對(duì)每一個(gè)節(jié)點(diǎn)都采用式(2)在其積分域上積分，可以得到所有離散節(jié)點(diǎn)的系統(tǒng)方程，將這些離散節(jié)點(diǎn)的系統(tǒng)方程組裝起來(lái)就能夠獲得問(wèn)題域的整體系統(tǒng)方程。

利用移動(dòng)最小二乘法(MLS)得到節(jié)點(diǎn)積分點(diǎn)支撐域內(nèi)的形函數(shù)，從而獲得問(wèn)題域的位移逼近函數(shù):

式中:ΦT(X)為根據(jù)MLS所得的形函數(shù)矩陣，uI為離散節(jié)點(diǎn)I的節(jié)點(diǎn)值，N為積分點(diǎn)的支撐域Ωs中的節(jié)點(diǎn)數(shù)。

根據(jù)彈性力學(xué)應(yīng)力－應(yīng)變關(guān)系有

式中:D為材料彈性矩陣，B為幾何矩陣，且

式中:ΓQt為積分域ΩQ和自然邊界的重合部分，N為自然邊界上的單位外法向矢量矩陣。

1.2剛度方程

把節(jié)點(diǎn)系統(tǒng)方程式(3)和式(4)代入節(jié)點(diǎn)控制方程式(5)，可以得到

式(6)可以簡(jiǎn)記為矩陣形式:

其中

這樣由每個(gè)離散節(jié)點(diǎn)的控制方程總裝成整個(gè)結(jié)構(gòu)系統(tǒng)的控制方程表達(dá)為

由于在MLS中，采用的節(jié)點(diǎn)支撐域都是緊支的，系統(tǒng)的總體剛度矩陣K是帶狀稀疏矩陣，可以減少計(jì)算工作量。根據(jù)文獻(xiàn)［13］對(duì)有限元法中罰因子的分析，建議罰因子α 105～12×E在范圍內(nèi)取值，E為楊氏彈性模量。

2　B樣條小波基

在MLPG法中，運(yùn)用最小二乘法(MLS)是求解位移場(chǎng)逼近函數(shù)的關(guān)鍵，從節(jié)點(diǎn)控制方程的推導(dǎo)過(guò)程來(lái)看，權(quán)函數(shù)的選定及其參數(shù)的確定對(duì)于計(jì)算的精度和結(jié)果的收斂至關(guān)重要。在無(wú)網(wǎng)格法中，一般要求權(quán)函數(shù)為緊支函數(shù)，在求解域內(nèi)具有連續(xù)可導(dǎo)、單調(diào)遞減和正態(tài)緊支等特性［14］。由于小波函數(shù)具有在不同分辨率水平上表達(dá)函數(shù)的變尺度能力，而其相應(yīng)的基函數(shù)(小波基)則在具有震蕩性、衰減性的同時(shí)還具有緊支性或近似緊支性［15－16］;另一方面B樣條函數(shù)具有分段光滑、局部可微和線性組合等特征［17］，所以B樣條小波函數(shù)的時(shí)頻局部特性最接近無(wú)網(wǎng)格法中常用作權(quán)函數(shù)的Gauss函數(shù)，并且它的緊支性優(yōu)于Gauss函數(shù)，這樣保證了在計(jì)算處理過(guò)程中的相位失真。所以本文將利用小波函數(shù)和B樣條函數(shù)的各自優(yōu)點(diǎn)，來(lái)構(gòu)造小波基函數(shù)，并以此構(gòu)造出來(lái)的B樣條小波基函數(shù)作為MLS法中的權(quán)函數(shù)。

依據(jù)式(9)的節(jié)點(diǎn)序列可以進(jìn)一步構(gòu)造在[0，1 ]區(qū)間內(nèi)的m階局部支撐嵌套的j尺度逼近空間V，其基函數(shù)則可以表達(dá)為

式中:Nm(x)為基數(shù)樣條函數(shù)。

局部域內(nèi)的小波緊支撐空間表示為

文章選用三次B樣條小波函數(shù)，即m=3。三次B樣條小波函數(shù)的基本表達(dá)式為［17］

3　本質(zhì)邊界條件

因無(wú)網(wǎng)格法形函數(shù)ΦT(X)不滿足克羅內(nèi)克條件，對(duì)于剛度方程式(7)中的節(jié)點(diǎn)I的參數(shù)u并非其真實(shí)位移，所以本質(zhì)邊界條件不能夠像有限元法那樣直接施加。Chen J S等提出了全轉(zhuǎn)換法(Full transformation method)［19］來(lái)處理其邊界條件。本文按照全轉(zhuǎn)換法的原理對(duì)剛度方程(控制方程)進(jìn)行修正，得到

這樣就可以得到節(jié)點(diǎn)的真實(shí)位移u－，然后可以直接施加本質(zhì)邊界條件。

4　算例

為了驗(yàn)證文章方法的正確性，本節(jié)給出船舶結(jié)構(gòu)中“船體梁”和“開(kāi)孔板”兩種典型結(jié)構(gòu)的計(jì)算分析。本文算例中材料的楊氏彈性模量為E=2.1× 105MPa，泊松比μ=0.3。

4.1船體梁

梁結(jié)構(gòu)是船體基本結(jié)構(gòu)之一。圖1為一端固支、一端自由的船體簡(jiǎn)化梁結(jié)構(gòu)。均布載荷p=10 N/m。

圖1　受均布載荷的船體梁Fig.1 Typical beam structure under uniform pressure

圖2(a)是根據(jù)本文無(wú)網(wǎng)格算法經(jīng)過(guò)4步自適應(yīng)生成的離散節(jié)點(diǎn)圖。初始離散為90個(gè)均布節(jié)點(diǎn)，積分背景網(wǎng)格為351個(gè)，經(jīng)過(guò)4步自適應(yīng)分析后，節(jié)點(diǎn)數(shù)為387個(gè)，背景網(wǎng)格為1 210個(gè)。圖2(b)為ANSYS生成的有限元網(wǎng)格單元。

圖2　無(wú)網(wǎng)格法與FEM節(jié)點(diǎn)與網(wǎng)格圖Fig.2 Arrangement of the node with Mesh-free and FEM

采用本文無(wú)網(wǎng)格方法，計(jì)算出算例梁的應(yīng)力云圖，如圖3(a)所示。通過(guò)與圖3(b)比較，應(yīng)力集中區(qū)域表現(xiàn)都是一致的。

圖3　無(wú)網(wǎng)格法與FEM計(jì)算應(yīng)力云圖Fig.3 The Mises-stress map with MLPG and FEM

圖4和圖5給出了在小波基尺度scale=2.5下按照本文無(wú)網(wǎng)格法和有限元ANSYS計(jì)算結(jié)果和解析解比較曲線(圖中FEM為有限元法解，ANS為解析解，MLPG為文章算法解，后面圖例一致)。為了簡(jiǎn)單起見(jiàn)，僅對(duì)比板Y方向中線截面上正應(yīng)力分布情況和板X方向中線截面上位移變化圖。從圖可以得出:1)FEM法和MLPG法的計(jì)算相對(duì)誤差范數(shù)都在5％以下;2)文章提出的方法比FEM方法更接近解析解?？梢?jiàn)本文方法具有很高的計(jì)算精度。

圖4　板X方向中線截面上位移變化比較圖Fig.4　Total deformation of the middle section along the X-direction

表1為在小波基不同尺度下的計(jì)算誤差值表(圖中Lu為位移相對(duì)誤差范數(shù)，Lσ為應(yīng)力相對(duì)誤差范數(shù)，后面相同)，根據(jù)表1可以繪制誤差曲線圖為圖6所示。從表1和圖6的變化趨勢(shì)可以得出結(jié)論:1)計(jì)算精度受尺度變化影響，也就是支撐域的大小影響到計(jì)算精度;2)進(jìn)一步提高文章方法的精度在于確定恰當(dāng)?shù)男〔ɑ叨?，可以由尺度三次B樣條函數(shù)的尺度函數(shù)來(lái)調(diào)節(jié)。

圖5　板X方向中線截面上正應(yīng)力變化圖Fig.5 Mises-stress of the middle section along the X-direction

表1　三次B樣條基函數(shù)MLPG法尺度－相對(duì)誤差范數(shù)表Table 1 The relative error norm with B-spline wavelet ％

圖6　不同尺度下相對(duì)誤差范數(shù)變化曲線圖Fig.6 The relative error norm with 3-order B-spline wavelet

4.2開(kāi)孔板

圖7為簡(jiǎn)化為兩端自由、另兩端承受均勻拉力的平直船體板，板中間開(kāi)有圓孔。此類型結(jié)構(gòu)在船舶結(jié)構(gòu)上也非常普遍，板側(cè)均布拉力q=1 MPa。

圖8(a)是根據(jù)本文無(wú)網(wǎng)格算法經(jīng)過(guò)四步自適應(yīng)生成的離散節(jié)點(diǎn)圖。初始離散為749個(gè)均布節(jié)點(diǎn)，積分背景網(wǎng)格為1 390個(gè)，經(jīng)過(guò)四步自適應(yīng)分析后，節(jié)點(diǎn)數(shù)為1 323個(gè)，背景網(wǎng)格為2 515個(gè)。圖8(b)為ANSYS生成的有限元網(wǎng)格單元。

圖7　中間開(kāi)孔平直板示意圖Fig.7 Flat plate with a circular opening hole

圖8　無(wú)網(wǎng)格法與FEM節(jié)點(diǎn)與網(wǎng)格圖Fig.8 Arrangement of the node with Mesh-free and FEM

采用本文無(wú)網(wǎng)格方法，計(jì)算出算例板的應(yīng)力云圖，如圖9(a)所示。通過(guò)與有限元軟件ANSYS計(jì)算的應(yīng)力云圖圖9(b)所示比較，應(yīng)力分布表現(xiàn)都是一致的。

圖9　無(wú)網(wǎng)格法與FEM計(jì)算應(yīng)力云圖Fig.9 The Mises-stress map with MLPG and FEM

進(jìn)一步比較本文無(wú)網(wǎng)格法與有限元ANSYS計(jì)算的結(jié)果的相似程度。為了簡(jiǎn)單起見(jiàn)，僅對(duì)比板Y方向中線截面上正應(yīng)力分布情況和板X方向中線截面上位移變化圖。結(jié)果如圖10和圖11所示。FEM法的位移相對(duì)誤差范數(shù)最大為2.5％、應(yīng)力誤差范數(shù)最大為2.4％，而文章方法的位移相對(duì)誤差范數(shù)為0.3％、應(yīng)力誤差范數(shù)最大為0.28％，可見(jiàn)文章方法的計(jì)算精度具有明顯的優(yōu)勢(shì)。

為了進(jìn)一步了解文章方法的正確性，下面給出在尺度為2.1沿板中線X方向節(jié)點(diǎn)使用FEM法和文章方法計(jì)算的位移相對(duì)誤差范數(shù)變化圖，如圖12所示。由圖和前面的分析可以得出結(jié)論:1)文章方法具有很好的計(jì)算精度;2)尺度選定對(duì)文章算法有很大的影響，不同尺度下的計(jì)算精度有很大差別(如算例1中圖6所示規(guī)律)。

圖10　板Y方向中線截面上正應(yīng)力變化圖Fig.10 Mises-stress of the middle section along the Y-direction

圖11　板X方向中線截面上位移變化圖Fig.11 Total deformation of the middle section along the X-direction

圖12　尺度2.1下的位移相對(duì)誤差范數(shù)變化曲線圖Fig.12 The relative error norm with 3-order B-spline wavelet and Gauss weight function

5　結(jié)論

文章提出的基于小波基的船舶結(jié)構(gòu)無(wú)網(wǎng)格分析技術(shù)由于擁有較好的自適應(yīng)性，在計(jì)算精度上比傳統(tǒng)的有限元法具有優(yōu)勢(shì)。采用有限元法(ANSYS)和提出的方法對(duì)船體梁和開(kāi)孔板進(jìn)行進(jìn)行實(shí)例計(jì)算，并對(duì)計(jì)算結(jié)果進(jìn)行比較和分析，可以得出結(jié)論:

1)文章方法具有很好的計(jì)算精度;

2)尺度選定對(duì)文章算法有很大的影響，不同尺度下的計(jì)算精度有很大差別;

3)當(dāng)3階B樣條小波函數(shù)作為權(quán)函數(shù)，結(jié)果權(quán)重函數(shù)比高斯權(quán)函數(shù)更好;

4)文章所提出的方法可以在小支撐域的情況下，達(dá)到良好的精度，即它可以實(shí)現(xiàn)更高的擬合。

參考文獻(xiàn):

［1］孫麗萍，李力波.船舶結(jié)構(gòu)有限元分析［M］.哈爾濱:哈爾濱工程大學(xué)出版社，2013:1－12.SUN Liping，LI Libo.Finite element analysis of ship structures［M］.Harbin:Harbin Engineering University Press，2013:1－12.

［2］BELYTSCHKO T，KRONGAUZ Y，ORGAN D，et al.Meshless methods:An overview and recent developments ［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1－4):3－47.

［3］LIU W K，HAO S，BELYTSCHKO T，et al.Multiple scale meshless methods for damage fracture and localization［J］.Computational Materials Science，1999，16(1/2/3/4):197－205.

［4］ATLURI S N，SHEN Shengping.The basis of meshless domain discretization:the meshless local Petrov－Galerkin(MLPG)method［J］.Advances in Computational Mathematics，2005，23(1/2):73－93.

［5］ODEN J T，DUARTE C A M，ZIENKIEWICZ O C.A new cloud－based hp finite element method［J］.Computer Methods in Applied Mechanics and Engineering，1998，153(1/2):117－126.

［6］LIU G R，GU Y T.Meshless local Petrov－Galerkin(MLPG)method in combination with finite element and boundary element approaches［J］.Computational Mechanics，2000，26(6):536－546.

［7］何沛祥，李子然，吳長(zhǎng)春.無(wú)網(wǎng)格與有限元的耦合在動(dòng)態(tài)斷裂研究中的應(yīng)用［J］.應(yīng)用力學(xué)學(xué)報(bào)，2006，23(2):195－198.HE Peixiang，LI Ziran，WU Changchun.Coupled finite element－element－free Galerkin method for dynamic fracture［J］.Chinese Journal of Applied Mechanics，2006，23(2):195－198.

［8］段念，王文珊，于怡青，等.基于FEM與SPH耦合算法的單顆磨粒切削玻璃的動(dòng)態(tài)過(guò)程仿真［J］.中國(guó)機(jī)械工程，2013，24(20):2716－2721.DUAN Nian，WANG Wenshan，YU Yiqing，et al.Dynamic simulation of single grain cutting of glass by coupling FEM and SPH［J］.China Mechanical Engineering，2013，24(20):2716－2721.

［9］JOHNSON G R，STRYK R A，BEISSEL S R，et al.An algorithm to automatically convert distorted finite elements into meshless particles during dynamic deformation［J］.International Journal of Impact Engineering，2002，27(10):997－1013.

［10］胡德安，韓旭，肖毅華，等.光滑粒子法及其與有限元耦合算法的研究進(jìn)展［J］.力學(xué)學(xué)報(bào)，2013，45(5):639－652.HU Dean，HAN Xu，XIAO Yihua，et al.Research developments of smoothed particle hydrodynamics method and its coupling with finite element method［J］.Chinese Journal of Theoretical and Applied Mechanics，2013，45(5):639－652.

［11］LIU G R.Meshfree methods:moving beyond the finite element method［M］.Boca Raton:CRC Press，2002:56－85.

［12］ATLURI S N，ZHU T.A new meshless local Petrov－Galerkin(MLPG)approach in computational mechanics［J］.Computational Mechanics，1998，22(2):117－127.

［13］ZIENKIEWICZ O C，TAYLOR R L.The finite element method［M］.4th ed.London:McGraw－Hill，1989:97－101.

［14］楊玉英，李晶.無(wú)網(wǎng)格Galerkin方法中權(quán)函數(shù)的研究［J］.塑性工程學(xué)報(bào)，2005，12(4):5－9.YANG Yuying，LI Jing.A study of weight function in element－free Galerkin method［J］.Journal of Plasticity Engineering，2005，12(4):5－9.

［15］LONG Shuyao，HU Dean.A study on the weight function of the moving least square approximation in the local boundary integral equation method［J］.Acta Mechanica Solida Sinica，2003，16(3):276－282.

［16］MALLAT S G.Multiresolution approximations and wavelet orthonormal bases of L2(R)［J］.Transactions of the American Mathematical Society，1989，315(1):69－87.

［17］秦榮.樣條無(wú)網(wǎng)格法［M］.北京:科學(xué)出版社，2012:25－46.QIN Rong.Spline meshless method［M］.Beijing:Science Press，2012:25－46.

［18］CHUI C K，QUAK E.Wavelets on a bounded interval ［M］//BRAESS D，SCHUMAKER L L.eds.Numerical Methods of Approximation Theory.Basel:Birkh?user Basel，1993，9:53－57.

［19］CHEN J S，PAN Chunhui，WU Chengtang，et al.Reproducing kernel particle methods for large deformation analysis of non－linear structures［J］.Computer Methods in Applied Mechanics and Engineering，1996，139(1/2/3/4):195－227.

Meshless analysis method of ship structures based on a B-spline wavelet

CHEN Jianping1，2，TANG Wenyong1，XU Manping2
(1.School of Naval Architecture，Ocean＆Civil Engineering，Shanghai Jiao Tong University，Shanghai 200240，China;2.School of Ship Engineering，Guangzhou Maritime Institute，Guangzhou 510725，China)

Abstract:To solve the problem of the high gradient adaptive analysis of the ship straight structure，a meshless local Petrov－Galerkin method based on a B－spline wavelet was proposed.The approximation function of the structural displacement field quantity was solved by employing the least squares method and the weighted residual method，and the governing equation and stiffness equation were established.Based on the meshless local Petrov－Galerkin method，an m－order B－spline function was used as the wavelet basis function to construct the approximation function of the ship structure displacement field，and a two－scale decomposition technology was used to decompose the high gradient component and the low scale component in the stress field.The high scale component was used to express the high gradient component in the stress field.Numerical examples show the solutions are good agreement between FEM－ANSYS and the proposed method，which verifies the validity of the presented method for analysis of the ship structures.

Keywords:ship structure analysis;approximation function;moving least square method;B－spline wavelet;scaling function;meshless local meshless patrov－galerkinmethod

通信作者:陳建平，E－mail:wchchenjp@ sina.com.

作者簡(jiǎn)介:陳建平(1973－)，男，副教授，博士.

基金項(xiàng)目:國(guó)家自然科學(xué)基金重點(diǎn)資助項(xiàng)目(51239007);中國(guó)博士后基金資助項(xiàng)目(2015M581622);廣東省自然科學(xué)基金資助項(xiàng)目(2014A030313792).

收稿日期:2014－09－05.網(wǎng)絡(luò)出版時(shí)間:2015－12－21.

中圖分類號(hào):U661.42

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):1006－7043(2016)01－0013－06

doi:10.11990/jheu.201409018

網(wǎng)絡(luò)出版地址:http://www.cnki.net/kcms/detail/23.1390.u.20151221.1613.044.html