基于統(tǒng)計(jì)方法的NOPD耗能機(jī)理定量分析

崔致遠(yuǎn)，吳九匯，陳花玲，李滌塵

(西安交通大學(xué) 機(jī)械工程學(xué)院，西安 710049)

非阻塞性微顆粒阻尼(NOPD)技術(shù)又稱粉體阻尼技術(shù)，是一種近些年發(fā)展起來(lái)的新型阻尼減振技術(shù)。它是在結(jié)構(gòu)振動(dòng)傳輸路徑上，加工一定數(shù)量的孔洞，在其中填充適當(dāng)數(shù)量的金屬或非金屬球形顆粒(顆粒直徑范圍為0.005～0.5 mm)。隨著結(jié)構(gòu)體的振動(dòng)，粉體顆粒之間不斷撞擊和摩擦，不但有動(dòng)量交換，而且能夠消耗系統(tǒng)的振動(dòng)能，達(dá)到減振的目的。NOPD具有很多優(yōu)點(diǎn)，比如:對(duì)原系統(tǒng)改動(dòng)很小、幾乎不產(chǎn)生附加質(zhì)量、顯著提高結(jié)構(gòu)的阻尼比、能適用于惡劣環(huán)境、減振性能不會(huì)隨時(shí)間降低等，具有廣泛的應(yīng)用前景［1－4］。

針對(duì) NOPD的研究，此前多處于實(shí)驗(yàn)分析階段［5－6］。隨著顆粒物質(zhì)力學(xué)的發(fā)展，散體單元法(DEM)也已開(kāi)始應(yīng)用于NOPD的減振機(jī)理研究;基于內(nèi)時(shí)理論的研究也有人涉及。散體單元法(DEM)是一種不連續(xù)數(shù)值方法模型，其優(yōu)點(diǎn)在于能夠考慮散體中實(shí)際顆粒的組成結(jié)構(gòu)，并能根據(jù)靜力學(xué)或動(dòng)力學(xué)原理研究單個(gè)顆粒及其總和的性質(zhì)，適用于模擬離散顆粒組合體在準(zhǔn)靜態(tài)或動(dòng)態(tài)條件下的變形過(guò)程［7－9］。文獻(xiàn)［10－11］結(jié)合NOPD的結(jié)構(gòu)特點(diǎn)，構(gòu)造了一種球體元模型，認(rèn)為NOPD的耗能機(jī)理分為2種:一種為沖擊耗能，另一種為摩擦耗能。影響沖擊耗能的主要因素為彈性恢復(fù)系數(shù)，影響摩擦耗能的主要因素為摩擦系數(shù)和法向作用力。文獻(xiàn)通過(guò)計(jì)算機(jī)仿真得到NOPD阻尼機(jī)理的一般性結(jié)論，初步分析了顆粒填充率、顆粒大小等因素對(duì)能量耗散的影響。文獻(xiàn)［12］也由DEM理論出發(fā)，應(yīng)用牛頓第二定律建立了單個(gè)顆粒的運(yùn)動(dòng)方程，提出了可由4個(gè)彈簧阻尼器和3個(gè)摩擦阻尼器代表的顆粒間的接觸簡(jiǎn)化模型，根據(jù)赫茲接觸理論分析顆粒與顆粒及顆粒與結(jié)構(gòu)之間的接觸關(guān)系，同時(shí)采用“盒式”計(jì)算的基本思想，簡(jiǎn)化了顆粒間接觸判別算法，該方法在計(jì)算顆粒數(shù)較大的問(wèn)題時(shí)，可節(jié)省大量計(jì)算時(shí)間。在上述理論指導(dǎo)下得出的計(jì)算機(jī)仿真結(jié)果顯示:系統(tǒng)響應(yīng)隨著質(zhì)量比的增加而減小，也隨著顆粒材料密度的增大而減小，與實(shí)驗(yàn)結(jié)果具有較好的一致性。但是，當(dāng)NOPD結(jié)構(gòu)的微顆粒數(shù)目較多(超過(guò)1萬(wàn)粒)，或者NOPD的結(jié)構(gòu)較為復(fù)雜時(shí)，該模型的計(jì)算效率及精度顯得不令人滿意。

內(nèi)蘊(yùn)時(shí)間理論(簡(jiǎn)稱內(nèi)時(shí)理論)由Valanis等［13］提出，并用于描述耗散材料粘塑性過(guò)程的理論，內(nèi)時(shí)理論是通過(guò)對(duì)由內(nèi)變量表征的材料內(nèi)部組織結(jié)構(gòu)不可逆變化所滿足的熱力學(xué)約束條件的研究，得到內(nèi)變量變化所必須滿足的規(guī)律，從而給出具體材料在具體條件下一條特定的不可逆熱力學(xué)變量的演變路徑。文獻(xiàn)［14］應(yīng)用內(nèi)時(shí)理論分析NOPD的結(jié)構(gòu)響應(yīng)，確定了材料的內(nèi)時(shí)本構(gòu)特性，建立起散粒體的增量型內(nèi)時(shí)本構(gòu)方程，并基于此給出了NOPD結(jié)構(gòu)的有限元?jiǎng)恿Ψ匠?，并用New mark方法對(duì)動(dòng)力方程進(jìn)行了數(shù)值計(jì)算，結(jié)論認(rèn)為NOPD阻尼對(duì)振幅較大結(jié)構(gòu)的薄弱模態(tài)有較好的減振效果，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果具有較好的一致性。應(yīng)用上述方法建立的NOPD模型，雖與實(shí)驗(yàn)數(shù)據(jù)有較好的一致性，也具有一定的工程指導(dǎo)意義，但上述研究并未得出NOPD能量耗散的定量規(guī)律，顆粒直徑、材料密度及顆粒流體積比(顆粒體積/孔洞體積)對(duì)能量耗散率的具體影響依然不清楚。

顆粒物質(zhì)理論中，將顆粒流分為彈性流(彈性準(zhǔn)靜態(tài)顆粒流)和慣性流(慣性碰撞顆粒流)，其研究對(duì)象為處于相對(duì)簡(jiǎn)單運(yùn)動(dòng)狀態(tài)的顆粒流。NOPD中顆粒的運(yùn)動(dòng)包含了擠壓、相對(duì)滑動(dòng)、碰撞等，顯然運(yùn)動(dòng)形式更為復(fù)雜，簡(jiǎn)單彈性流及慣性流的知識(shí)并不能直接應(yīng)用到NOPD。本文借鑒了湍流的統(tǒng)計(jì)處理方法，從研究顆粒之間的能量耗散出發(fā)，以NOPD顆粒群整體作為研究對(duì)象，建立了顆粒流運(yùn)動(dòng)方程及能量方程，得到了NOPD單位質(zhì)量、單位時(shí)間的能量耗散率與顆粒大小及顆粒流體積比等參數(shù)的相互關(guān)系，所得理論結(jié)果與已有實(shí)驗(yàn)數(shù)據(jù)具有良好的一致性。

1 NOPD振動(dòng)模型的統(tǒng)計(jì)處理

本節(jié)基于湍流的耗能統(tǒng)計(jì)模型，根據(jù)顆粒流類(lèi)流體性質(zhì)，從顆粒流一般本構(gòu)關(guān)系出發(fā)，借鑒局部各向同性湍流的耗能模型，得到了NOPD能量耗散率及能譜密度的表達(dá)式。

流體中，當(dāng)雷諾數(shù)超過(guò)某臨界值時(shí)，層流變得不穩(wěn)定，并開(kāi)始向湍流過(guò)渡。湍流場(chǎng)中，湍動(dòng)量及湍動(dòng)能在雷諾應(yīng)力作用下由均流傳遞到大渦，再由大渦逐級(jí)傳遞到小渦，最終在小渦中由于粘性作用耗散為熱。Kolmogorov［15］由此不僅提出局部各向同性的概念，還提出速度場(chǎng)結(jié)構(gòu)函數(shù)D(r)的概念來(lái)描述該類(lèi)湍流速度起伏強(qiáng)度的統(tǒng)計(jì)特征，并提出兩個(gè)著名假設(shè):① 湍流能量傳遞過(guò)程中，能量耗散率ε與運(yùn)動(dòng)粘性系數(shù)ν是決定能量傳遞的兩個(gè)特征量，并通過(guò)無(wú)因次分析可確定湍流場(chǎng)的特征長(zhǎng)度η和特征速度V;② 在慣性副區(qū)(r?η，r為湍流場(chǎng)中任意兩點(diǎn)間的距離，η為湍流特征尺寸)，湍流場(chǎng)的縱向速度關(guān)聯(lián)函數(shù)與運(yùn)動(dòng)粘性系數(shù)ν無(wú)關(guān)。

處于靜止或低速運(yùn)動(dòng)狀態(tài)下的顆粒流，各顆粒間相對(duì)運(yùn)動(dòng)不明顯，可認(rèn)為各點(diǎn)處的速度近似等于一平均速度，表現(xiàn)出與一般流體相似的性質(zhì)(類(lèi)流體性)。而NOPD處于工作狀態(tài)時(shí)，其顆粒流在外界振動(dòng)激勵(lì)下，各顆粒間發(fā)生明顯的相對(duì)運(yùn)動(dòng)，顆粒流內(nèi)部發(fā)生對(duì)流［16］，隨激振頻率的增高，對(duì)流現(xiàn)象越明顯，不同速度的顆粒之間發(fā)生相互擠壓、摩擦、碰撞，部分顆粒將不再跟隨顆粒流整體一起運(yùn)動(dòng)，顆粒的運(yùn)動(dòng)狀態(tài)較運(yùn)動(dòng)初期顯得混亂無(wú)序，當(dāng)NOPD達(dá)到工作平衡狀態(tài)時(shí)，以近似平均速度運(yùn)動(dòng)的顆粒流將被具有不同速度的顆粒取代，大量顆粒在振動(dòng)激勵(lì)下以各自速度往復(fù)運(yùn)動(dòng)，呈現(xiàn)出一定的周期性。顆粒對(duì)流運(yùn)動(dòng)實(shí)際上是具有相同運(yùn)動(dòng)速度的顆粒組成的“相關(guān)層”之間有規(guī)則的剪切運(yùn)動(dòng)［17］，顆粒流內(nèi)的動(dòng)量傳遞便是由剪切運(yùn)動(dòng)導(dǎo)致的速度波動(dòng)引起［18］。Taguchi［19－20］曾用顆粒湍流模型，用數(shù)值模擬方法研究湍流，其研究結(jié)果顯示:當(dāng)振動(dòng)強(qiáng)度超過(guò)某個(gè)臨界值時(shí)，顆粒發(fā)生流化，出現(xiàn)與完全發(fā)展的湍流類(lèi)似的渦結(jié)構(gòu)，并且，顆粒流中的能譜密度存在與湍流中類(lèi)似的規(guī)律(即kolmogorov用于描述局部各向同性湍流的－5/3定律)。由此，振動(dòng)顆粒流可看作局部各向同性。所以，本文中將采用湍流局部各向同性模型描述振動(dòng)顆粒流。

顆粒流動(dòng)的一般本構(gòu)關(guān)系由文獻(xiàn)［23］用連續(xù)介質(zhì)力學(xué)中一般方法得到:

其中，C為顆粒間粘滯力，φ0為顆粒材料靜止內(nèi)摩擦角，p為接觸壓力，d為顆粒粒徑，ρs為顆粒密度，D為變形速度張量，I為張量D的不變量，kτ，kp，為與顆粒體積比ψ、顆粒流發(fā)生剪切變形的臨界體積比以及顆粒彈性恢復(fù)系數(shù)有關(guān)的應(yīng)力系數(shù)，J2’=1/2trD2，kτ1=kp1tanφ1，kτ2=kp2tanφ2，φ1，φ2分別定義為滑動(dòng)摩擦角和碰撞摩擦角。NOPD顆粒間存在法向粘滯力及切向摩擦力，但法向粘滯力的影響較小，一般不予考慮(C=0)，由于顆粒對(duì)流實(shí)際上是有規(guī)則的簡(jiǎn)單剪切運(yùn)動(dòng)，即為簡(jiǎn)單剪切流V1=U(y)，V2=V3=0，本構(gòu)方程可簡(jiǎn)化為:

其中速度梯度(dU/dy)零次項(xiàng)由顆粒之間靜態(tài)支撐作用引起;速度梯度線性項(xiàng)由顆粒之間的相對(duì)滑動(dòng)和擠壓作用引起，該項(xiàng)可被看作耗散項(xiàng)，顆粒間通過(guò)相互擠壓、摩擦耗散外界能量，因此與能量耗散有關(guān)的參數(shù)為速度梯度線性項(xiàng)系數(shù)kp1ρsd3/2g1/2與kτ1ρsd3/2g1/2，其中kp1，kτ1為無(wú)因次量;速度梯度二次項(xiàng)是由顆粒間的碰撞和擴(kuò)散作用引起，該項(xiàng)為對(duì)流擴(kuò)散項(xiàng)，其作用是使能量在顆粒間重新分配，并不改變能量總和。

令ν=d3/2g1/2，并引入表征顆粒流特征參數(shù):特征速度u與特征尺度η，通過(guò)無(wú)因次分析得到由ε與ν表示的特征參數(shù)表達(dá)式:

當(dāng)兩點(diǎn)速度與其連線同向時(shí)，用Bdd(r)表示該兩點(diǎn)的速度關(guān)聯(lián)函數(shù)，稱為縱向關(guān)聯(lián)函數(shù)。依據(jù)相似假設(shè)，在局部各向同性區(qū)域內(nèi)，Bdd(r)與等效運(yùn)動(dòng)粘性系數(shù) ν 無(wú)關(guān)［15，24］，得到:

其中，βdd為(r/η)的普適函數(shù)，C1為常數(shù)，由實(shí)驗(yàn)確定。

能譜密度也有類(lèi)似定律，在振動(dòng)強(qiáng)度較高的情況下有:

由式(4)、式(5)可發(fā)現(xiàn)，NOPD能量耗散率及能譜密度的表達(dá)式形式與湍流理論中Kolmogorov假設(shè)的相關(guān)結(jié)論一致，將通過(guò)實(shí)驗(yàn)數(shù)據(jù)分析驗(yàn)證該結(jié)論。

2 NOPD顆粒流的能量關(guān)系

從N－S方程出發(fā)，結(jié)合顆粒流一般本構(gòu)關(guān)系，建立適應(yīng)顆粒流動(dòng)量守恒方程(廣義N－S方程)及能量守恒方程，分析NOPD工作狀態(tài)下的能量轉(zhuǎn)移、分布及耗散等狀態(tài)間的關(guān)系。

確定偏應(yīng)力時(shí)，顆粒流一般本構(gòu)關(guān)系中的速度梯度零次項(xiàng)是不重要的［25］，由無(wú)粘條件，可將本構(gòu)關(guān)系式(1)寫(xiě)為:

將式(6)代入N－S方程，令:

可得到不可壓縮條件下顆粒流的動(dòng)量守恒方程(廣義N－S方程)和連續(xù)方程:

考慮連續(xù)方程，整理(9)式右邊第一項(xiàng)，得:

其中Ⅰ項(xiàng)為脈動(dòng)運(yùn)動(dòng)中粘性切力做的功，即脈動(dòng)速度產(chǎn)生的動(dòng)能平均輸送項(xiàng)，表示脈動(dòng)對(duì)平均流的影響;Ⅱ項(xiàng)為脈動(dòng)運(yùn)動(dòng)中的能量耗損項(xiàng)，用ε表示。

3 NOPD顆粒流的能量關(guān)系數(shù)據(jù)分析

根據(jù)以上文中得到的NOPD能量耗散率表達(dá)式，結(jié)合已有實(shí)驗(yàn)數(shù)據(jù)，畫(huà)出能量耗散率ε隨顆粒流體積比及顆粒大小等因素的變化曲線，將所得能量耗散率的變化規(guī)律與已有的相關(guān)結(jié)果進(jìn)行對(duì)比。處于湍動(dòng)平衡狀態(tài)的顆粒流，雖然整體上并不處于簡(jiǎn)單剪切流狀態(tài)，但可認(rèn)為其中某小部分顆粒群處于簡(jiǎn)單剪切運(yùn)動(dòng)狀態(tài)。由無(wú)粘條件及簡(jiǎn)單剪切流假設(shè)，公式(10)中Ⅱ項(xiàng)(能量耗散率ε)可表示如下:

圖1 三種不同顆粒直徑(d)的顆粒流在du/dy=100 s－1時(shí)，能量耗散率 ε隨顆粒體積比ψv的變化曲線Fig.1 Mean energy dissipation(ε)curve(log － log)for different values of volume ratio(ψv)，when d=1.8mm，d=1.1mm，d=0.5mm;du/dy=100 s－1.

圖2 三種不同體積比(ψv)顆粒流在du/dy=100 s－1時(shí)，能量耗散率ε隨顆粒直徑d的變化曲線Fig.2 Mean energy dissipation(ε)curve for different values of particle diameter(d)，when ψv=0.56，ψv=0.54，ψv=0.49，du/dy=100 s－1.

由圖1可見(jiàn)，在相同顆粒大小、相同速度梯度情況下，顆粒流體積比越高顆粒間發(fā)生的碰撞及摩擦越多，從而使NOPD的能量耗散率隨體積比的增加而平滑上漲;當(dāng)顆粒流體積比超過(guò)60%時(shí)，NOPD的能量耗散率迅速增加，并趨于其最大值。在研究顆粒流體積比對(duì)NOPD阻尼效果影響的文獻(xiàn)中，也有相同的變化規(guī)律，即振動(dòng)的耗散時(shí)間在顆粒流體積比超過(guò)60%時(shí)急劇減少，并趨于其最小值［27］。圖 2為在相同體積比時(shí)，NOPD的能量耗散率隨顆粒尺寸的增大而增加，體積比較高時(shí)，顆粒尺寸的增加會(huì)導(dǎo)致更加顯著的能量耗散率的變化。此前的研究結(jié)果也證實(shí)了該變化規(guī)律的正確性［6］。上述對(duì)比證明了利用湍流耗能模型分析NOPD耗能問(wèn)題的可行性和正確性。

圖3所示，能譜密度E與頻率f的變化成反比。當(dāng)振動(dòng)強(qiáng)度較高時(shí)，顆粒間才會(huì)發(fā)生對(duì)流現(xiàn)象，此時(shí)顆粒相關(guān)層中的剪切運(yùn)動(dòng)比較低振動(dòng)強(qiáng)度時(shí)更劇烈，從而導(dǎo)致更多的能量耗散。該現(xiàn)象與湍流由慣性區(qū)向耗能區(qū)的轉(zhuǎn)變過(guò)程極為相似。

圖3 ψv=0.56，d=1.1 mm 的干顆粒流，能譜密度E隨頻率f的變化曲線Fig.3 Energy spectral density(E)curve for different values of frequency(f)，when d=1.1mm，ψv=0.56.

4 結(jié)論

本文提出了一個(gè)分析NOPD能量耗散機(jī)理的定量模型，并得到了滿意的理論結(jié)果，結(jié)論表明:顆粒直徑越大、顆粒流體積比越高，NOPD的能量耗散率就越高;相同能量耗散率情況下，能譜密度隨著波數(shù)的增加而降低。能量耗散率的分析結(jié)果與參考文獻(xiàn)中的實(shí)驗(yàn)數(shù)據(jù)有較好的一致性，說(shuō)明湍流耗能模型對(duì)NOPD的耗能過(guò)程具有適用性。

該理論模型能更好的適應(yīng)大量的顆粒及復(fù)雜結(jié)構(gòu)條件下NOPD的耗能分析，對(duì)工程應(yīng)用有一定的指導(dǎo)作用。該模型的建立，為更精確地分析NOPD耗能機(jī)理以及更好地指導(dǎo)工程應(yīng)用開(kāi)辟了一條新的思路。

［1］Xu Z W，Wang M Y，Chen T N.A particle damper for vibration and noise reduction［J］.Journal of Sound and Vibration，2004，207(4 －5):1033 －1040.

［2］Panossian H V.Structural damping enhancement via nonobstructive particle damping technique［J］. Journalof Vibration and Acoustics，1992，114:101 －105.

［3］Panossian H V，Kovac B，Rackl R.Composite honeycomb treatment via non-obstructive particle damping(NOPD)［R］.Proceedings ofthe45thAIAA/ASME/ASCE/AHS/ASC Structures;Structural Dynamics ＆ Materials Conference，2004，1689.

［4］Panossian H V.Ehrgott R.Non-obstructive particle damping(NOPD)treatment optimization for composite honeycomb panels［R］.Proceedings of the 48th AIAA/ASME/ASCE/AHS/ASC Structures:Structural Dynamics and Materials Conference，2007，2047.

［5］Xu Z W，Wang M Y，Chen T N.Particle damping for passive vibration suppression:numericalmodeling and experimentalinvestigation［J］. JournalofSound and Vibration，2005，279:1097 －1120.

［6］Xu Z W，Wang M Y，Chen T N.An experimental study of particle damping for beams and plates［J］.Journal of Vibration and Acoustics，2004，126:141 －148.

［7］Mao K M，Wang M Y，Xu Z W，et al.DEM simulation of particle damping［J］.Powder Technology，2004，142:154－165.

［8］Wong C X，Daniel M C，Rongong J A.Energy dissipation prediction of particle dampers［J］.Journal of Sound and Vibration，2009，319(1 －2):91 －118.

［9］Saeki M.Impact damping with granular materials in a horizontally vibrating system［J］.Journal of Sound and Vibration，2002，251(1):153－161.

［10］MaoK ，Wang M，Ding H，et al.Simulation and characterization of particle damping in transient vibrations［J］.Journal of Vibration and Acoustics，2004，126(2):202－211.

［11］Chen T N，Mao K ，Huang X Q，et al.Dissipation mechanism of non-obstructive particle damping using discrete element method［J］.In:Proceedings of SPIE International Symposium on Smart Structures and Materials.Damping and Isolation，2001，4331:294 －301.

［12］Hu L，Huang Q，Liu Z.A non-obstructive particle damping model of DEM［J］.International Journal of Mechanics and Materials Design，2008，4:45 －51.

［13］ Valanis K C，F(xiàn)an J.A numerical algorithm for endochronic plasticity and comparison with experiment［J］.Computers ＆Structures，1984，19(5):717 －724.

［14］王 煒，黃協(xié)清，陳天寧，等.內(nèi)蘊(yùn)時(shí)間理論用于NOPD結(jié)構(gòu)相應(yīng)計(jì)算的研究［J］.力學(xué)學(xué)報(bào)，2003，35(2):246－252.

［15］ Kolmogorov A N.The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers［J］.Proceedings:Mathematical and Physical Sciences，Turbulence and Stochastic Process:Kolmogorov＇s Ideas 50 Years On，1991，434(1890):9 －13.

［16］ Eshuis P，Meer D V D，Alam M，et al.Onset of convection in strongly shaken granular matter［J］.Physical Review Letters，2010，104(3):038001.

［17］ McNamara S，Luding S.Energy flows in vibrated granular media［J］.Physical Review E，1998，58:813 －822.

［18］Gioia G，Ott-Monsivais E，Hill K M.Fluctuating velocity and momentum transfer in dense granular flows［J］.Physical Review Letters，2006，96(13):138001.

［19］ Taguchi Y H.Turbulent flow in vibrated bed of powder:new target to investigate turbulent flow［J］.Fractals，1993，1(4):1080－1085.

［20］ Taguchi Y H.K-5/3 Power spectrum in Powder-turbulent flow in a vibrated bed:numerical results［J］.Europhysics Letters，1993，24(3):203 －209.

［21］ Taguchi Y H.Powder turbulence:direct onset of turbulent flow［J］.Journal de Physique II 2，1992:2103 －2114.

［22］ Losert W，Bocquet L，Lubensky T C，et al.Particle dynamics in sheared granular matter［J］.Physical Review Letters，2000，85:1428 －1431.

［23］王光謙，熊 剛，方紅衛(wèi).顆粒流動(dòng)的一般本構(gòu)關(guān)系［J］.中國(guó)科學(xué)E輯:技術(shù)科學(xué)，1998，28(3):282－288.

［24］Kolmogorov A N.Dissipation of energy in the locally isotropic turbulence［J］.Proceedings:Mathematical and Physical Sciences，Turbulence and Stochastic Process:Kolmogorov＇s Ideas 50 Years On，1991，434(1890):15 －17.

［25］韓式方.非牛頓流體連續(xù)介質(zhì)力學(xué)［M］.成都:四川科學(xué)技術(shù)出版社，1988.

［26］ 陸厚根.粉體工程導(dǎo)論［M］.上海:同濟(jì)大學(xué)出版社，1993.

［27］Wakasawa Y，Hashimoto M，Marui E.The damping capacity improvement of machine tool structures by balls packing［J］.International Journal of Machine Tools＆ Manufacture，2004，44(14):1527－1536.