泵用不可壓縮流體密封剛度系數(shù)分析

張盟，王曉放，徐勝利，萬(wàn)學(xué)麗

(1.大連理工大學(xué) 能源與動(dòng)力學(xué)院，遼寧大連116023;2.大連深藍(lán)泵業(yè)有限公司，遼寧 大連116031)

泵葉輪入口輪蓋與進(jìn)口導(dǎo)流管間存在動(dòng)靜間隙。部分高壓流體經(jīng)葉輪出口間隙外泄，并重新回流至泵入口，這股回流既消耗主泵的功率，也干擾主流場(chǎng)流動(dòng)，同時(shí)減小有效通流面積，降低泵的流動(dòng)效率和性能。為了盡量降低這種泄漏，在動(dòng)靜間隙上設(shè)置非接觸式不可壓縮流體密封。

非接觸式密封技術(shù)廣泛應(yīng)用在泵等旋轉(zhuǎn)機(jī)械中，它能夠有效控制旋轉(zhuǎn)部件與靜止部件間的泄漏。對(duì)于轉(zhuǎn)子動(dòng)力系統(tǒng)，密封會(huì)提供附加的剛度和阻尼，這對(duì)轉(zhuǎn)子的穩(wěn)定性造成一定的影響［1］。

Childs等［2-4］對(duì)控制體方法不斷改進(jìn)。Arghir等［5-7］發(fā)展了 CFD方法，并計(jì)算動(dòng)力特性系數(shù)。Benckert等［8］做了大量的關(guān)于動(dòng)力特性的實(shí)驗(yàn)，并證明密封的交叉剛度是由密封周向流動(dòng)引起的，他們?cè)趯?shí)驗(yàn)中測(cè)量了不同類型密封的直接剛度，發(fā)現(xiàn)較長(zhǎng)密封的直接剛度為負(fù)值。Leong等［9］蒸汽輪機(jī)迷宮密封做了大量試驗(yàn)，結(jié)果與Benckert等的測(cè)量結(jié)果很吻合，多數(shù)密封直接剛度為負(fù)，少量短密封為正。Mihai等［10］發(fā)現(xiàn)進(jìn)出口壓差較小時(shí)，出口有回堵現(xiàn)象，動(dòng)力特性系數(shù)中的直接剛度系數(shù)出現(xiàn)負(fù)值，影響轉(zhuǎn)子的對(duì)中效應(yīng)，轉(zhuǎn)子穩(wěn)定性差。

國(guó)內(nèi)何立東等［11-13］用實(shí)驗(yàn)和數(shù)值方法研究動(dòng)力特性系數(shù)。模擬仿真方面，孫婷梅等［14-15］利用CFD有限元軟件Fluent計(jì)算迷宮密封三維流場(chǎng)，研究了偏心率、入口預(yù)旋、渦動(dòng)速度對(duì)密封動(dòng)力特性的影響，密封直接剛度維持在負(fù)值范圍，他們計(jì)算出密封周向壓力分布曲線，但是沒(méi)有分析產(chǎn)生負(fù)直接剛度的原因。

本文應(yīng)用數(shù)值模擬結(jié)合工程實(shí)際，利用CFD有限元軟件Fluent計(jì)算LNG泵的節(jié)流襯套和口環(huán)密封的動(dòng)力特性系數(shù)，模擬在LNG泵不同轉(zhuǎn)速工況下，不可壓縮流體密封的動(dòng)特性變化，通過(guò)分析比較密封圓周上壓力和速度分布，研究等截面環(huán)形密封產(chǎn)生負(fù)直接剛度的原因和影響因素。

1 密封模型

1.1 渦動(dòng)轉(zhuǎn)子數(shù)學(xué)模型

本文假設(shè)整體計(jì)算域?yàn)橥牧?，湍流模型采用?biāo)準(zhǔn)的k-ε模型［16］，近壁面采用標(biāo)準(zhǔn)壁面函數(shù):

式中:Gk表示由平均速度梯度產(chǎn)生的湍動(dòng)能項(xiàng)，Gb表示由浮力產(chǎn)生的湍動(dòng)能項(xiàng)，YM表示由可壓縮湍流中，耗散率的波動(dòng)項(xiàng)，方程常數(shù)項(xiàng)C1ε=1.44，C2ε=1.92，Cμ=0.09，還有 σk=1.0，σε=1.3 分別是湍動(dòng)能k和湍動(dòng)能耗散率ε的湍流普朗特?cái)?shù)，Sk和Sε是自定義源項(xiàng)，標(biāo)準(zhǔn)k-ε湍流模型是半經(jīng)驗(yàn)公式。

1.2 渦動(dòng)轉(zhuǎn)子動(dòng)力學(xué)模型

本文采用旋轉(zhuǎn)坐標(biāo)系，在轉(zhuǎn)子中心定義坐標(biāo)系，轉(zhuǎn)子與靜子的相對(duì)位置不變，坐標(biāo)系變換把非定常問(wèn)題轉(zhuǎn)換成定常問(wèn)題。圖1所示偏心轉(zhuǎn)子在靜子中渦動(dòng)受力，e為轉(zhuǎn)子偏心距，本文假設(shè)轉(zhuǎn)子繞靜子中心以圓形軌跡渦動(dòng)，渦動(dòng)半徑是轉(zhuǎn)子偏心距e，渦動(dòng)角速度是Ω，轉(zhuǎn)子的旋轉(zhuǎn)角速度是ω，轉(zhuǎn)子的旋轉(zhuǎn)角速度與渦動(dòng)角速度比值定義為渦動(dòng)比，渦動(dòng)軌跡的徑向力和切向力為Fr、Ft。

圖1 偏心轉(zhuǎn)子在靜子中渦動(dòng)受力示意圖Fig.1 Force on eccentric rotor in stator

1.3 動(dòng)力學(xué)方程

當(dāng)轉(zhuǎn)子受到小擾動(dòng)，以小圓形軌跡繞靜子中心渦動(dòng)，轉(zhuǎn)子受到的水動(dòng)力可以由剛度、阻尼和慣性系數(shù)的反對(duì)稱矩陣和轉(zhuǎn)子的位移、速度、加速度的線性關(guān)系表示，矩陣中包括 6個(gè)獨(dú)立參數(shù)(K，k，C，c，M，m)，K為直接剛度，k為交叉剛度，C為主阻尼，c為交叉阻尼，M為主慣性系數(shù)，m為耦合慣性系數(shù)。

在不同的渦動(dòng)比下，積分轉(zhuǎn)子表面壓力生成轉(zhuǎn)子渦動(dòng)軌跡的徑向力和切向力Fr、Ft，徑向力和切向力的二階方程:

本文求解密封線性動(dòng)力特性，偏心距選取徑向間隙的10%，Moore等［17］驗(yàn)證過(guò)此值在模擬小軌跡渦動(dòng)的準(zhǔn)確性。

為了求解特性系數(shù)，至少選取3個(gè)渦動(dòng)比，為了提高計(jì)算精度，一般選取6個(gè)渦動(dòng)比，對(duì)結(jié)果進(jìn)行線性回歸計(jì)算。

1.4 計(jì)算域和邊界條件

圖1(b)所示轉(zhuǎn)子和靜子相對(duì)位置的軸向視圖，在計(jì)算域入口指定總壓，出口指定靜壓。旋轉(zhuǎn)坐標(biāo)系下，轉(zhuǎn)子壁面相對(duì)角速度為ω-Ω，密封壁面相對(duì)角速度為-Ω。給定湍流強(qiáng)度和水力直徑。在轉(zhuǎn)子和密封壁面為標(biāo)準(zhǔn)壁面條件，采用速度無(wú)滑移條件，流動(dòng)絕熱。

1.5 求解方法

求解器選擇分離的隱式求解器，采用有限體積法離散控制方程。連續(xù)方程、動(dòng)量方程和能量方程的離散格式為二階迎風(fēng)格式，湍動(dòng)能和耗散率方程采用一階迎風(fēng)格式。

1.6 網(wǎng)格

在偏心狀態(tài)下生成三維結(jié)構(gòu)化六面體網(wǎng)格，偏心距選取間隙的10%。圖2為槽道式密封間隙和槽道二維網(wǎng)格。為了提高求解精度，在近壁面增加了網(wǎng)格密度，相鄰兩節(jié)點(diǎn)間距離之比為1.1。Hirano等［7］驗(yàn)證了密封空腔和間隙網(wǎng)格密度獨(dú)立性，證明了本文網(wǎng)格密度計(jì)算結(jié)果的可靠性。

圖2 槽道式密封間隙和槽道二維網(wǎng)格Fig.2 Mesh of land and groove for grooved seal

2 計(jì)算結(jié)果與分析

2.1 等截面環(huán)形密封數(shù)值模擬

圖1(a)為等截面環(huán)形密封示意圖，對(duì)不同間隙等截面環(huán)形密封的不同轉(zhuǎn)速工況進(jìn)行數(shù)值模擬，轉(zhuǎn)速nr=500～10 000 r/min，額定轉(zhuǎn)速n=2 980 r/min，模擬介質(zhì)為液態(tài)LNG，其他參數(shù)為:內(nèi)徑Di=53.2 mm，外徑Do=53.5 mm，密封長(zhǎng)度L=50 mm，L/Di=0.94(長(zhǎng)密封)，介質(zhì)進(jìn)出口壓差 ΔP=0.005 MPa，環(huán)境壓力P0=0.1MPa，動(dòng)力粘度 μ=0.000 125 kg/(m·s)。

計(jì)算1 000～10 000 r/min轉(zhuǎn)速工況下，1倍間隙和2倍間隙(磨損)等截面環(huán)形密封剛度特性系數(shù)，分析轉(zhuǎn)速對(duì)剛度特性系數(shù)的影響，如圖3所示，等截面環(huán)形密封直接剛度為負(fù)值，1倍間隙直接剛度的絕對(duì)值大于2倍間隙的絕對(duì)值，且直接剛度的絕對(duì)值隨著轉(zhuǎn)速的增大而增大，1倍間隙的增大幅度比2倍間隙大。1倍間隙的交叉剛度大于2倍間隙的交叉剛度。轉(zhuǎn)速增大時(shí)，交叉剛度增大，且1倍交叉剛度的增幅大于2倍交叉剛度。

如圖4所示，在5 000r/min工況下，等截面環(huán)形密封最大間隙和最小間隙的速度和壓力特性曲線，等截面環(huán)形密封介質(zhì)壓力沿軸向下降，斜率和速度有關(guān)。密封周向間隙不均勻，介質(zhì)受高速旋轉(zhuǎn)的轉(zhuǎn)子的摩擦力影響，在最小間隙處形成較最大間隙更大的周向速度，因此在最小間隙處形成較最大間隙處更大的速度。如圖4(a)在密封進(jìn)口處，最大間隙和最小間隙的速度激增，且最小間隙的加速度大于最大間隙，這部分動(dòng)能增量消耗更多壓力能，如圖4(b)密封進(jìn)口處，最小間隙壓力下降幅度明顯大于最大間隙;在密封中部，速度基本保持定值，最大間隙和最小間隙壓降斜率一致;在接近密封出口，最大間隙處壓力急劇下降，最小間隙處緩慢下降，同時(shí)下降到出口壓力，在此過(guò)程最大間隙和最小間隙的速度相應(yīng)幅度的增長(zhǎng)。在密封軸向長(zhǎng)度上，最小間隙處壓力始終小于最大間隙處，因此表明直接剛度為負(fù)值，且直接剛度系數(shù)與壓力特性曲線包含的面積有關(guān)。

圖3 等截面環(huán)形密封剛度特性曲線Fig.3 Stiffness characteristic of annular seal

如圖5(a)所示，在極坐標(biāo)系下，在2 950 r/min工況下，1倍間隙和2倍間隙等截面環(huán)形密封2/5軸長(zhǎng)處，圓周上各角度速度，介質(zhì)在密封環(huán)內(nèi)成螺旋軌跡流動(dòng)，在圓周方向上最小間隙受到最大流動(dòng)阻力，產(chǎn)生節(jié)流效應(yīng)，1倍間隙密封在270°位置(最小間隙)速度達(dá)到最大，如圖5(b)在174°位置壓力最大，在334°壓力最小，因此以 174°、334°為分界點(diǎn)，從174°到334°在轉(zhuǎn)速方向上提供較大的壓力勢(shì)，促進(jìn)介質(zhì)加速流動(dòng);然而從334°到174°轉(zhuǎn)速方向上提供反向壓力勢(shì)阻止介質(zhì)流動(dòng)，在88°位置速度最低。在同樣的壓差邊界條件下，2倍間隙泄漏大于1倍間隙，因此各角度速度大于1倍間隙。

圖4 等截面環(huán)形密封(5 000 r/min)最大間隙和最小間隙的速度和壓力特性曲線Fig.4 Velocity and pressure of annular seal maximum and minimum clearance in the condition of 5 000 r/min

圖5(b)中，壓力特性橢圓中心向上偏離，表征密封產(chǎn)生負(fù)直接剛度;向左偏離，表征產(chǎn)生正交叉剛度;偏離極坐標(biāo)系中心越大，表征剛度值越大。1倍間隙密封在較窄的流道內(nèi)受到較大的阻力，所以壓力特性橢圓中心偏心較大，1倍間隙的交叉剛度大于2倍間隙，直接剛度絕對(duì)值稍大于2倍間隙。因此，密封從正常間隙磨損成2倍間隙可以減小負(fù)直接剛度的不穩(wěn)定影響。

圖6所示等截面環(huán)形密封不同轉(zhuǎn)速圓周壓力分布對(duì)比圖，如圖6(a)，隨轉(zhuǎn)速增大，壓力特性橢圓中心向左上方偏離距離逐漸增大，表征負(fù)直接剛度和交叉剛度的絕對(duì)值逐漸增大，穩(wěn)定性下降。

如圖6(b)，不同轉(zhuǎn)速下壓力特性點(diǎn)對(duì)比圖，特性點(diǎn)包含的面積與直接剛度的絕對(duì)值成正比，壓力特性點(diǎn)面積的長(zhǎng)軸線斜率的絕對(duì)值與交叉剛度的絕對(duì)值成正比。隨著轉(zhuǎn)速增大，直接剛度的絕對(duì)值逐漸增大，交叉剛度逐漸增大。10 000 r/min工況下負(fù)直接剛度絕對(duì)值最大，1 000 r/min工況下負(fù)直接剛度絕對(duì)值最小;10 000 r/min工況下正交叉剛度最大，1 000 r/min工況下正交叉剛度最小。

如圖6(c)，1 000 r/min和500 r/min壓力特性點(diǎn)對(duì)比圖，1 000 r/min壓力特性點(diǎn)出現(xiàn)交叉點(diǎn)，500 r/min壓力特性交叉點(diǎn)向z軸正方向移動(dòng)，直接剛度由負(fù)值變?yōu)檎?。因此，推測(cè)隨著轉(zhuǎn)速的繼續(xù)降低，直接剛度會(huì)逐漸增大。推論高轉(zhuǎn)速是造成負(fù)剛度現(xiàn)象的因素之一。

圖5 等截面環(huán)形密封(2 950 r/min)1倍間隙和2倍間隙的速度和壓力特性橢圓Fig.5 Velocity and pressure circumferential distribution of annular seal normal clearance and twice clearance in the condition of 2 950 r/min

圖6 等截面環(huán)形密封不同轉(zhuǎn)速圓周壓力分布對(duì)比圖Fig.6 Comparison of pressure distribution for annular seal over a range of rotating speed

2.2 槽道式密封數(shù)值模擬

圖7為槽道式密封示意圖，對(duì)槽道式密封的不同轉(zhuǎn)速工況進(jìn)行數(shù)值模擬，轉(zhuǎn)速nr=1 000～7 000 r/min，額定轉(zhuǎn)速n=2 980 r/min，模擬介質(zhì)為液態(tài)LNG，其他參數(shù)為:內(nèi)徑Di=194.65 mm，外徑Do=195.05 mm，密封長(zhǎng)度L=28.5 mm，槽道數(shù)量N=6，槽道深度h=2 mm，槽道寬度w=2 mm，介質(zhì)進(jìn)出口壓差 ΔP=0.37 MPa，環(huán)境壓力P0=0.1 MPa，動(dòng)力粘度 μ=0.000 125 kg/m·s。

計(jì)算在1 000～10 000r/min轉(zhuǎn)速工況下，1倍間隙和2倍間隙(磨損)槽道式密封剛度特性系數(shù)，分析轉(zhuǎn)速對(duì)剛度特性系數(shù)的影響，如圖8所示，槽道式密封1倍間隙直接剛度小于2倍間隙直接剛度，且直接剛度隨著轉(zhuǎn)速呈現(xiàn)緩慢下降的趨勢(shì)。1倍間隙的交叉剛度大于2倍間隙的交叉剛度。轉(zhuǎn)速增大時(shí)，交叉剛度增大，且1倍交叉剛度的增幅大于2倍交叉剛度。

圖9所示槽道式密封圓周壓力特性點(diǎn)對(duì)比圖，由于槽道式密封有較大密封腔室，所以密封圓周壓力和速度分布比等截面環(huán)形密封更均勻，最大間隙和最小間隙之間壓差較等截面環(huán)形密封更小，更不會(huì)出現(xiàn)最小間隙壓力比最大間隙小的現(xiàn)象。在相同壓差邊界條件下，2倍間隙槽道式密封的泄漏量較1倍間隙大，速度也較大，因此壓力較1倍間隙小。對(duì)于槽道式密封圓周平均壓力越小，壓力特性點(diǎn)包含的面積越大，槽道式密封直接剛度越大。2倍間隙槽道式密封的直接剛度大于1倍間隙，交叉剛度小于1倍間隙。

圖10所示為2倍間隙槽道式密封不同轉(zhuǎn)速下壓力特性點(diǎn)對(duì)比圖。對(duì)于較大間隙槽道式密封，隨著轉(zhuǎn)速增大，圓周平均壓力值逐漸增大，且趨于均勻，壓力特性點(diǎn)包含的面積逐漸縮小，直接剛度逐漸減小，交叉剛度逐漸增大，但增大幅度較小。1 000 r/min工況下直接剛度最大，7 000 r/min工況下直接剛度最小。7 000 r/min工況下交叉剛度最大，1 000 r/min工況下交叉剛度最小。

圖7 槽道式密封示意圖Fig.7 Grooved seal

圖8 槽道式密封剛度特性曲線Fig.8 Stiffness characteristic of grooved seal

圖9 槽道式密封(7 000 r/min)圓周壓力特性點(diǎn)對(duì)比圖Fig.9 Pressure change with Z-coordinate for grooved seal at 7 000 r/min

圖10 槽道式密封不同轉(zhuǎn)速圓周壓力特性點(diǎn)對(duì)比圖Fig.10 Pressure change with Z-coordinate for grooved seal over a range of rotating speed

3 結(jié)論

1)等截面環(huán)形密封和槽道式密封的1倍間隙直接剛度小于2倍間隙，1倍間隙交叉剛度大于2倍間隙，密封經(jīng)過(guò)磨損泄漏量增大，穩(wěn)定性卻提高。

2)轉(zhuǎn)速使等截面環(huán)形密封和槽道式密封的直接剛度下降，交叉剛度增大，造成穩(wěn)定性下降，且隨轉(zhuǎn)速不穩(wěn)定隱患成倍增長(zhǎng)。

3)對(duì)于不可壓縮流體，較長(zhǎng)的等截面環(huán)形密封容易產(chǎn)生負(fù)直接剛度，槽道式密封直接剛度通常為正值，比等截面環(huán)形密封更加穩(wěn)定，工程上建議采用槽道式密封控制不可壓縮流體泄漏。

4)對(duì)于等截面環(huán)形密封，高轉(zhuǎn)速容易造成負(fù)剛度現(xiàn)象的重要因素，但對(duì)于旋轉(zhuǎn)機(jī)械，為避免這種現(xiàn)象，工程上建議增大間隙以減小轉(zhuǎn)子不穩(wěn)定隱患。

影響負(fù)剛度的因素很多，在結(jié)構(gòu)參數(shù)上，密封的L/D比值是產(chǎn)生負(fù)剛度的重要的結(jié)構(gòu)參數(shù)，本文計(jì)算的等截面環(huán)形密封L/Di=0.94屬于長(zhǎng)密封，容易產(chǎn)生負(fù)剛度且與短密封動(dòng)特性有很大差別。針對(duì)長(zhǎng)密封的結(jié)構(gòu)參數(shù)對(duì)動(dòng)特性影響還有待于進(jìn)一步研究。

［1］BLACK H F，COCHRANE E A.Leakage and hybrid bearing properties of serrated seals in centrifugal pumps［C］//The 6th International Conference on Fluid Sealing.Munich，Germany，1973:G5-61-G5-70.

［2］CHILDS D W，SCHARRER J K.An Iwatsubo based solution for labyrinth seals:a comparison to experimental results［J］.ASME J Eng Gas Turbines Power，1986，108:325–331.

［3］IWATSUBO T，SHENG B.Evaluation of dynamic characteristics of parallel grooved seals by theory and experiment［C］//Proceedings of the Third IFToMM International Conference on Rotordynamics.Lyon，F(xiàn)rance，1990:313-318.

［4］FLORJANCIC S.Annular seals of high energy centrifugal pumps:a new theory and full scale measurement of rotordynamic coefficients and hydraulic friction factors［D］.Zurich:Swiss Federal Institute of Technology，1990:26-63.

［5］ARGHIR M，F(xiàn)RENE J.Rotordynamic coefficients of circumferentially-grooved liquid seals using the averaged Navier-Stokes equations［J］.ASME J Tribology，1997，119:556-567.

［6］MOORE J J.Three-dimensional CFD rotordynamic analysis of gas labyrinth seal［J］.ASME Journal of Vibrations and Acoustics，2003，125:427-433.

［7］HIRANO T，GUR Z，KIRK R G.Application of computational fluid dynamics analysis for rotating machinery-part II:labyrinth seal analysis［J］.ASME J Eng Gas Turbines Power，2005，127:820-826.

［8］BENEKERT H，WACHTER J.Flow induced spring coefficients of labyrinth seals for applications in turbomachinery［C］//Rotordynamic Instability Problems in High-Performance Turbomachinery.Texas:NASA，1980:189-212.

［9］LEONG Y，BROWN R D.Circumferential pressure distributions in a model labyrinth seal［C］//Rotordynamic Instability Problems in High-Performance Turbomachinery.Texas:NASA，1982:223-241.

［10］MIHAI A，CYRIL D，JEAN F.The Lomakin effect in annular gas seals under choked flow conditions［J］.ASME Journal of Tribology，2007，129(4):1028-1034 .

［11］何立東，高金吉，金琰，等.三維轉(zhuǎn)子密封系統(tǒng)氣流激振的研究［J］.機(jī)械工程學(xué)報(bào)，2003，39(3):100-104.HE Lidong，GAO Jinji，JIN Yan，et al.Study on gas flow-induced vibration for a three-dimensional rotor-seal system［J］.Chinese Journal of Mechanical Engineering，2003，39(3):100-104.

［12］王煒哲，劉應(yīng)征，葉春，等.迷宮密封-轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)特性的試驗(yàn)測(cè)量和數(shù)值模擬［J］.機(jī)械工程學(xué)報(bào)，2007，43(3):22-27.WANG Weizhe，LIU Yingzheng，YE Chun，et al.Experimental measurement and numerical simulation of dynamics of labyrinth seal—rotor combination［J］.Chinese Journal of Mechanical Engineering，2007，43(3):22-27.

［13］曹浩，楊建剛，郭瑞，等.密封動(dòng)力特性系數(shù)試驗(yàn)識(shí)別方法及影響因素分析［J］.機(jī)械工程學(xué)報(bào)，2011，47(9):85-89.CAO Hao，YANG Jiangang，GUO Rui，et al.Experimental identification method and influence factor analysis of seal dynamic characteristic［J］.Chinese Journal of Mechanical Engineering，2011，47(9):85-89.

［14］孫婷梅，鄭水英.基于FLUENT迷宮密封動(dòng)力特性分［J］.流體機(jī)械，2008，136(108):24-28.SUN Tingmei，ZHENG Shuiying.Dynamic coefficient analysis of labyrinth seal based on FLUENT［J］.Fluid Machinery，2008，136(108):24-28.

［15］孫婷梅.迷宮密封流場(chǎng)及其動(dòng)力特性計(jì)算［D］.杭州:浙江大學(xué)，2008:30-48.SUN Tingmei.Flow and dynamic coefficients of labyrinth seals［D］.Hangzhou:Zhejiang University，2008.

［16］LAUNDER B E，SPALDING D B.Lectures in mathematical models of turbulence［M］.London:Academic Press，1972:32-126.

［17］MOORE J J，PALAZZOLO A B.CFD comparison to 3D laser anemometer and rotordynamic force measurements for grooved liquid annular seals［J］.Journal of Tribology-transactions of The ASME，1999，121(2):307-314.