張鵬舉++吳曉冬++李博++吳恭興

摘要： 建立500 m水深的采油樹的流場簡化模型，對采油樹流場在給定的2種工況下進行分析，評定流體流動和分布特點并計算流阻.分析結(jié)果表明：采油樹本體在2種工況下的壓力和速度變化較均勻，只是在變直徑處或直角彎管處會出現(xiàn)湍流現(xiàn)象.

關鍵詞： 采油樹； 流場； 流阻； 湍流； 有限元； FLUENT

中圖分類號： TE53；TB115.1文獻標志碼： B

0引言

21世紀是海洋的世紀，隨著海洋的開發(fā)與利用，海洋石油成為熱點話題.我國南海海域?qū)掃_2×106 km2，是世界上四大海洋油氣聚集中心之一.我國海洋石油儲量約為2.30×1010～3.00×1010 t，天然氣儲量3.38×1014 m3.[12]目前，國際海洋石油工程界普遍認可的深水定義是水深300 m.[3]我國水深300 m以上的海域達1.537×106 km2[4]，但大部分尚未被勘探，因此，深海石油的勘探將成為我國未來海洋資源開發(fā)的主要方向.

作為海洋油氣開發(fā)不可或缺的設備，水下采油樹在水下生產(chǎn)系統(tǒng)中起著至關重要的作用.[58]但是，由于我國深水石油的開采起步晚，技術不全面，深水采油樹的開發(fā)受到限制.長期以來，水下采油樹的關鍵技術被美國FMC和Cameron以及挪威Aker Kvaerner Subsea 三大廠商所壟斷，這3家企業(yè)占有世界采油樹市場90%以上的份額.由于水下采油樹的結(jié)構復雜，對材料性能和密封技術的要求很高，控制系統(tǒng)和閥等單元部件容易出現(xiàn)問題.目前，國內(nèi)海洋石油的裝備主要靠進口，因此對水下采油樹相關技術進行國產(chǎn)化研究，對突破國外技術封鎖有重大意義.[911]

1CAE分析簡介

1.1分析目的和要求

水下采油樹可分為立式和臥式2種，兩者的主要區(qū)別在于油管懸掛器和閘閥閥組的安裝位置不同.[12]本文研究的水下采油樹為臥式采油樹，其特點是油管掛在采油樹本體內(nèi)，采油樹的閥組位于油管掛的側(cè)面，生產(chǎn)油管和油管掛的安裝要后于采油樹自身.

采用通用CFD軟件FLUENT 6.3對采油樹流場進行CAE分析，評定流體流動和分布特點，并計算流阻.CAE分析和結(jié)果要滿足合適的API 6A的溫度要求：上游18～121 ℃，下游46～121 ℃.

1.2采油樹工況

采油樹以節(jié)流閥作為整棵樹的上下游分界，以流進節(jié)流閥的一端為上游，流出節(jié)流閥的一端為下游.上、下游油管流道由130.175 0 mm（518英寸）變?yōu)?03.187 5 mm（4116英寸）.

在采油樹實際工作過程中，流道內(nèi)一般為氣液混合相（CO2含量較高，基本沒有H2S），其設計使用壽命為20 a，設計水深為500 m.在流體特別是氣體流過節(jié)流閥時，由于壓差的存在，將產(chǎn)生絕熱節(jié)現(xiàn)象，即焦耳湯普生效應，節(jié)流閥及節(jié)流閥下游的設備會受此低溫效應的影響，因此上游主要設備（指與流體接觸的設備）額定設計溫度為18～121 ℃，下游設備額定設計溫度為46～121 ℃.

采油樹的額定壓力為69.0 MPa，測試壓力為103.5 MPa.實際上，節(jié)流閥的節(jié)流作用及流道孔徑由大變小導致下游壓力變小，也就是說，在上游負荷達到額定工作壓力的時候，下游設備并沒有達到額定工作壓力.

1.3分析內(nèi)容和假設

根據(jù)給定的工作環(huán)境，分析各組件關鍵部位在各種工況下流場的變化規(guī)律.

以當前的計算機計算速度和內(nèi)存，直接求解NS方程非常困難，故采用雷諾平均方程.

1.4計算方法

可壓縮流體流動和傳熱的控制方程可用雷諾平均NS方程表示，把流動變量放入連續(xù)性方程和動量方程，并且取一段時間的平均值，即

ρt+xi（ρui）=0 （1）

t（ρui）+xj（ρuiuj）=-pxi+

xjμuixj+ujxi-23δijulxl+xj（-ρu′iu′j）（2）

此方程適用于變密度的氣液混合流動，其中-ρu′iu′j稱為雷諾壓力，可利用Boussinesq假設把雷諾壓力與平均速度梯度連續(xù)起來，即-ρu′iu′j=μiuiuj+ujxi-23ρk+μtuixiδij（3）Boussinesq假設用在SA模型，kε模型和kω模型中，采用標準的kε湍流方程描述流體在管道內(nèi)的流動狀態(tài)，該方程主要基于湍流動能k和擴散率ε，即

t（ρk）+xi（ρkixi）=

xjμ+μiσkkxj+Gk+Gb-ρε-YM+Sk（4）

t（ρε）+xi（ρεμi）=xjμ+μiσεεxj+

C1εεk（Gk+C3εGk）-C2ερε2k+Sε （5）

式中：ρ為流體密度，kg/m3；μ為流體動力黏度；C1ε，C2ε和C3ε為經(jīng)驗常數(shù)，取值分別為1.44，1.92和0.99；σk和σε分別為k和ε的普朗特數(shù)，取值分別為1.0和1.3.[13]

此模型包含低雷諾數(shù)影響和可壓縮影響，適用于混合邊界層計算和受壁面限制的流動計算.

1.5采油樹整體分析

采油樹主要組成部分見圖1和2.H4連接器用于采油樹本體和井口的連接.油管掛坐于采油樹本體，并與采油樹本體的生產(chǎn)主閥和環(huán)空主閥相通.油管掛垂直中心孔以中心堵頭密封.油管掛往上為內(nèi)置樹帽（仍位于采油樹本體內(nèi)），再往上是垃圾帽，防止落入泥沙等.采油樹本體內(nèi)的環(huán)空主閥連接環(huán)空翼閥閥組，環(huán)空翼閥閥組連接環(huán)空跨接管.環(huán)空跨接管和采油樹本體的生產(chǎn)主閥都接入生產(chǎn)翼閥閥組.通過連接管1，生產(chǎn)翼閥閥組和雙孔連接器流入口連接，流入口經(jīng)過變徑〔130.175 0 mm（518英寸）變?yōu)?03.187 5 mm（4116英寸）〕進入RPM模塊內(nèi)的節(jié)流閥，然后又流回雙孔連接器的出口，經(jīng)過生產(chǎn)隔離閥，通過連接管2將直徑101.6 mm（4英寸）管道連接器連接起來.endprint

圖 1采油樹總裝圖

Fig.1Christmas tree assembly diagram

圖 2采油樹總裝圖剖面

Fig.2Profile of christmas tree of assembly diagram

在采油樹內(nèi)，流體的流通路線見圖3，通過H4連接器將采油樹本體和井口鎖緊固定，通過井下的油管將流體引入油管掛，再通過采油樹上的各類閥、閥組、連接管和連接器最終通過直徑101.6 mm（4英寸）管道連接器連接到外接管線.

圖 3流體流通路線

Fig.3Fluid circulation route

2CAE流場分析

2.1流體域模型

采油樹流體域模型見圖4.抽取流道壁并進行適當簡化，如忽略螺紋、小孔、小插拴和倒角等.圖 4采油樹流體域模型

Fig.4Fluid domain model of christmas tree

2.2網(wǎng)格劃分

網(wǎng)格劃分見圖5，共劃分2 381 521個四面體網(wǎng)格.由Gambit網(wǎng)格質(zhì)量檢查功能知：最差網(wǎng)格歪斜度為0.784 991，小于標準值0.97，所以所劃分的網(wǎng)格質(zhì)量較好，所有網(wǎng)格在整個流體計算域內(nèi)分布均勻、整齊，所以，用此網(wǎng)格進行計算能夠滿足工程精度要求，具體的網(wǎng)格質(zhì)量檢查見圖5c.

a）采油樹網(wǎng)格劃分

b1）節(jié)流閥b2）油管掛

b）節(jié)流閥和油管掛本體局部網(wǎng)格放大

c）采油樹網(wǎng)格質(zhì)量檢查

圖 5采油樹網(wǎng)格劃分和質(zhì)量檢查

Fig.5Christmas tree mesh and quality check

2.3邊界條件設定

計算域入口采用質(zhì)量流量入口，出口采用壓力出口，其他位置設置為固壁條件.認為管道中流動介質(zhì)為完全甲烷氣體.壓力差分格式采用標準離散差分格式，動量方程、動能方程和湍動能耗散率均采用2階迎風差分格式，通過SIMPLE算法耦合求解速度和壓力方程.

2.4計算工況

環(huán)境溫度：取海底深度為-340 m以下，中部深度約-170 m，平均環(huán)境溫度見表1.

表 1環(huán)境溫度

Tab.1Temperature of environment℃項目空氣海水海面中部海底最高36.031.423.519.5最低15.922.017.09.5

根據(jù)已知資料得到采油樹的最大日產(chǎn)量為6×105 m3，因此得到甲烷氣體在標準大氣壓下15.6 ℃時的體積流量Q=6×10524×3 600=6.944 m3/s （6）質(zhì)量流量M=ρ標Q （7）式中：M為質(zhì)量流量，kg/s；ρ標為甲烷氣體在標準大氣壓15.6 ℃時的密度，ρ標=0.677 6 kg/m3.具體計算過程見式（3）.計算工況見表2.

表 2計算工況

Tab.2Calculation conditions工況質(zhì)量流量

入口/（kg/s）出口

壓力/MPa入口

溫度/°C出口

溫度/°C環(huán)境

溫度/°C工況14.706 3.81032019工況24.70612.11032019

2.5設計壓力

采油樹額定工作壓力為69.0 MPa，測試壓力為103.5 MPa.

2.6分析計算

1）選擇FLUENT求解器.選擇三維單精度求解器進行分析.

2）網(wǎng)格的相關操作.讀入Gambit軟件生成的網(wǎng)格文件；檢查網(wǎng)格，確認最小網(wǎng)格體積小于零，否則需重新劃分網(wǎng)格；設置計算區(qū)域大小、模型和實際幾何尺寸的單位換算.

3）選擇計算模型.定義求解器，指定計算模型：選擇能量方程，選擇流態(tài)為湍流，湍流模式為kε兩方程模型，設置模型參數(shù).

4）設置操作環(huán)境.

5）定義流體的物理性質(zhì).將流入采油樹中的流體看作純甲烷氣體，將材料的物理性質(zhì)從數(shù)據(jù)庫中調(diào)出，并給定2種工況下的密度.

6）設置邊界條件.設置流體計算區(qū)域，設置質(zhì)量流量入口邊界條件和壓力出口邊界條件，其他邊界條件設置為固壁，設置入口溫度和出口溫度（見表2）.

7）求解方法的設置.設置求解參數(shù)，進行初始化；打開殘差圖，設置精度為0.000 01，在迭代計算時可動態(tài)顯示殘差，保存計算文件；設置迭代500次，進行迭代計算.

8）計算結(jié)果顯示.對給定的2種工況分別進行計算，完成后提取壓力云圖和速度云圖，見圖6～9，計算結(jié)果匯總見表3.

a）壓力云圖

b）節(jié)流閥和計量閥處壓力云圖局部放大

圖 6工況1采油樹壓力云圖

Fig.6Pressure contours of christmas tree under working condition 1

由圖6可知：最大壓力Pmax=24.1 MPa，最小壓力Pmin=3.2 MPa，壓差ΔP=20.9 MPa；由于節(jié)流閥和計量閥的作用，整個采油樹的壓力出現(xiàn)3段明顯的變化.由圖7可知：最大速度umax=739 m/s，最小速度umin=0；流體在經(jīng)過節(jié)流閥過程中速度增大，最大速度出現(xiàn)在節(jié)流閥和計量閥處，其他位置速度都很低.由圖8可知：最大壓力Pmax=24.4 MPa，最小壓力Pmin=7.95 MPa，壓差ΔP=16.45 MPa；由于節(jié)流閥的節(jié)流作用，壓力發(fā)生急劇變化.由圖9可知：最大速度umax=494 m/s，最小速度umin=0；由于節(jié)流閥的節(jié)流作用，流體在經(jīng)過節(jié)流閥時速度急劇增加.

a）速度云圖endprint

b）速度云圖局部放大

圖 7工況1采油樹速度云圖

Fig.7Velocity contours of christmas tree under working condition 1

a）壓力云圖

b）壓力云圖局部放大

圖 8工況2采油樹壓力云圖

Fig.8Pressure contours of christmas tree under working condition 2

a）速度云圖

b）速度云圖局部放大

圖 9工況2采油樹速度云圖

Fig.9Velocity contours of christmas tree under working condition 2

表 3計算結(jié)果

Tab.3Calculation results工況最大壓力/

MPa最小壓力/

MPa最大速度/

（m/s）最小速度/

（m/s）工況13.846 73.790 770.722 50工況212.138 712.093 121.594 50

3流阻計算

計算甲烷氣體在采油樹內(nèi)流動過程中密度隨壓力和溫度的變化.

由表3可知，甲烷氣體在流動過程中壓力變化不大，溫度變化也不大，所以密度基本不變.只需算在2種工況下的甲烷氣體密度，即可算出壓力損失.

真實氣體狀態(tài)方程為P=R（T+273.15）V-b-aV2 （8）式中：P為壓力，MPa；V為摩爾體積，m3/kmol；R為通用氣體常數(shù)或摩爾氣體常量，各種氣體的通用氣體常數(shù)相同，R=8.314 J/（mol·°C）；T為攝氏溫度；a為反應分子相互作用強度的常數(shù)，與物質(zhì)的性質(zhì)有關，一般由試驗獲得，對于甲烷氣體，a=0.228 589 m3/kmol；b為與分子體積有關的常數(shù)，與物質(zhì)的性質(zhì)有關，一般由試驗獲得，對于甲烷氣體b=0.042 71.

甲烷氣體在壓力12.1和3.8 MPa，溫度103 ℃時的摩爾體積可以根據(jù)式（8）進行計算，因此密度可以根據(jù)式（9）獲得.ρ甲烷=MmolV （9）式中：Mmol為甲烷氣體在相應壓力和溫度下的摩爾質(zhì)量.計算12.1和3.8 MPa時甲烷氣體的密度ρ12.1=66.147 kg/m3，ρ3.8=20.793 kg/m3.

實際流體在運動過程中必然會產(chǎn)生能量損失，在伯努利方程中將單位流體的機械能損失用hw表示，稱為水頭損失.根據(jù)產(chǎn)生水頭損失的流動邊界條件和機理不同，水頭損失可分為沿程水頭損失和局部水頭損失.

流體在進行均勻流動時，由于黏性作用，過流斷面上的流速分布不均勻，相鄰兩流層間存在相對運動，從而使流體流層之間及流體與邊界之間存在切應力，形成流動阻力.這種在均勻流段上產(chǎn)生的流動阻力稱為沿程阻力，由于沿程阻力做功而引起的水頭損失稱為沿程水頭損失，其計算公式為hf=fLQm/db （10）式中：f為沿程阻力因數(shù)，與管道的種類和材質(zhì)有關，一般由試驗獲得，此處根據(jù)相關資料取5 650；L為管道長度，m；Q為體積流量，m3/s，根據(jù)式（6）進行計算；m為流量指數(shù)，與管道的類型有關，一般由試驗獲得，此處根據(jù)相關資料取1.85；d為管徑，m；b為管徑指數(shù)，與管道的類型有關，一般由試驗獲得，此處根據(jù)相關資料取5.04.[14]

流體在流動的局部區(qū)域，如流體流經(jīng)管道的彎頭、突擴、突縮和閘門等處，由于固體邊界的急劇改變而引起速度分布的變化，甚至使主流脫離邊界，形成漩渦區(qū)從而產(chǎn)生的阻力稱為局部阻力.由于局部阻力做功引起的水頭損失稱為局部水頭損失，其計算公式為 hj=ξu22g （11）式中：ξ為局部阻力因數(shù)，雷諾數(shù)對ξ的影響很小，可以忽略不計，一般情況下ξ只取決于流道的邊界形式變化，可以根據(jù)經(jīng)驗選取；u為損失后的速度，m/s；g為重力加速度，m/s2.

由甲烷氣體在2種工況下的密度，可以根據(jù)質(zhì)量流量和管徑計算流速， M=ρuS （12）式中：M為質(zhì)量流量；ρ為甲烷氣體在相應壓力和溫度下的密度，kg/m3；u為流速，m/s；S為管道橫截面積，m2.

以直徑101.6 mm（4英寸）彎管1為例說明流阻計算過程.流體介質(zhì)完全為甲烷氣體，壓力為3.8 MPa，溫度為103.15 ℃.

將彎管1從油氣進口至油氣出口劃分為小段，依次編號為34，35，36和37，見圖10.34段為直管，只需計算沿程損失.體積流量Q=0.226 3 m3/s，直徑d34=0.102 m，管長L34=0.390 m，代入公式得沿程損失hf，34=0.040 m.37段為直管，同理Q=0.226 3 m3/s，d37=0.102 m，管長L37=0.652 m，代入公式得沿程損失hf，37=0.067 m.35和36段都是彎管，局部損失為主要部分，而沿程阻力相對局部阻力而言可以忽略.質(zhì)量流量M=4.706 kg/s，直徑d35=d36=0.101 m，彎管截面積S35=S36=0.008 2 m2.在壓力為3.8 MPa，溫度為103.15 ℃時，速度u=27.435 m/s.局部阻力因數(shù)根據(jù)參考資料分別取ξ35=0.30和ξ36=0.18.代入公式得局部損失hj，35=11.509 m，hj，36=6.905 m.

圖 10直徑101.6 mm（4英寸）彎管1

Fig.10Pipe 1 of 101.6 mm（4 m） diameter

根據(jù)上例，可以得到直徑101.6 mm（4英寸）彎管1的總水頭損失，匯總見表4.

表 4壓力為3.8 MPa，溫度為103.15 ℃時彎管1總損失

Tab.4Total pressure loss of pipe 1 under 3.8 MPa pressure and 103.15 ℃ temperature劃分段34353637總損失直徑d/mm102.0102.0截面積S/m20.008 20.008 2長度L/m0.390.39流量Q/（m3/s）0.226 30.226 3速度v/m/s27.43527.435沿程阻力因數(shù)5 6505 650局部阻力因數(shù)0.300.18沿程水頭損失/m0.0400.0670.107局部水頭損失/m11.5096.90518.414endprint

同理，采油樹整體的壓力損失計算結(jié)果見表5和6，采油樹內(nèi)流體流動過程見圖4.其中壓力損失ΔP=ρgh （12）式中：h為總水頭損失，單位m.由表6可知，當管道中流動的介質(zhì)完全是甲烷氣體時候，最大的壓力損失為1.596 3 MPa.

表 5工況1采油樹總流阻計算結(jié)果

Tab.5Total flow resistance calculation results of christmas tree under condition 1構件沿程水頭損失/m局部水頭損失/m總水頭頭損失/m壓力損失/MPa油管掛本體0.056 7.7267.7820.002直徑127.0 mm（5英寸）彎管0.112 36.48136.5930.007雙孔連接器0.6806 370.5166 371.1961.298直徑101.6 mm（4英寸）彎管10.12418.41418.5380.004生產(chǎn)隔離閥0.04300.0430直徑101.6 mm（4英寸）管20.41519.81520.230.004總損失1.4306 452.9526 454.3821.315

表 6工況2采油樹總流阻計算結(jié)果

Tab.6Total flow resistance calculation results of christmas tree under condition 2構件沿程水頭損失/m局部水頭損失/m總水頭頭損失/m壓力損失/MPa油管掛本體0.0061.0891.0950.000 7127.0 mm（5英寸）彎管0.0135.1205.1330.003 2雙孔連接器0.0792 450.2372 450.3161.588 4直徑101.6 mm（4英寸）彎管10.0142.7882.8020.001 8生產(chǎn)隔離閥0.00500.0050直徑101.6 mm（4英寸）管20.0483.3993.4470.002 2總損失0.1652 462.6332 462.7981.596 3

4結(jié)束語

利用FLUENT對采油樹內(nèi)流體流動特點和壓力的變化進行模擬，并計算流阻.結(jié)果表明：2種工況的壓力和速度變化較均勻，只是在變直徑處或直角彎管處會出現(xiàn)湍流現(xiàn)象.由流阻計算可知壓力損失約1.596 3 MPa.參考文獻：

[1]楊致政， 張瀛， 胡啟國. 近海邊際油田勘探開發(fā)中若干關鍵問題探討[J]. 特種油氣藏， 2005， 12（5）： 59.

YANG Zhizheng， ZHANG Ying， HU Qiguo. Approach to key points in exploration and development of offshore marginal oilfields[J]. Special Oil & Gas Reservoirs， 2005， 12（5）： 59.

[2]郭西水， 張林， 肖良平， 等. 深水完井技術研究[J]. 重慶科技學院學報（自然科學版）， 2011， 13（4）： 7476.

GUO Xishui， ZHANG Lin， XIAO Liangping， et al. Research on deep water well completion technology[J]. J Chongqing Univ Sci & Technol （Nat Sci）， 2011， 13（4）： 7476.

[3]馬彥金， 荊穎， 王亮. 對我國未來幾種深海油氣田開發(fā)新方案的探究[C]// 左其華， 竇希萍. 第十五屆中國海洋（岸）工程學術討論會論文集（上）. 北京： 海洋出版社， 2011： 522524.

[4]方華燦. 對我國深海油田開發(fā)工程中幾個問題的淺見[J]. 中國海洋平臺， 2006， 21（4）： 18.

FANG Huacan. A few superficial views about some problems of engineering of oil eield development in deep water in china[J]. China Offshore Platform， 2006， 21（4）： 18.

[5]任克忍， 王定亞， 周天明， 等. 海洋石油水下裝備現(xiàn)狀及發(fā)展趨勢[J]. 石油機械， 2008， 36（9）： 151153.

REN Keren， WANG Dingya， ZHOU Tianming， et al. Status and development of offshore oil and underwater equipment[J]. China Petroleum Machinery， 2008， 36（9）： 151153.

[6]楊進， 劉書杰， 周建良， 等. 深水石油鉆采工程模擬試驗裝置的研制[J]. 石油機械， 2011， 39（8）： 13.

YANG Jin， LIU Shujie， ZHOU Jianliang， et al. development of the simulation experimental device for deepwater petroleum drilling and production engineering[J]. China Petroleum Machinery， 2011， 39（8）： 13.

[7]趙宏林， 程華榮， 田紅平， 等. 深水采油樹油管懸掛器下放工具關鍵參數(shù)設計[J]. 石油機械， 2014， 42（4）： 1619.

ZHAO Honglin， CHENG Huarong， TIAN Hongping， et al. Design of the key parameters for the running tool of subsea christmas tree tubing hanger[J]. China Petroleum Machinery， 2014， 42（4）： 1619.endprint

[8]龔銘煊， 劉再生， 段夢蘭， 等. 深海水下采油樹安放過程分析與研究[J]. 石油機械， 2013， 41（3）： 5054.

GONG Mingxuan， LIU Zaisheng， DUAN Menglan， et al. Research on the running installment process of deep sea underwater christmas tree[J]. Petroleum Machinery， 2013， 41（3）： 5054.

[9]劉文霄， 鄧平， 金弢.海洋水下立式采油樹用油管懸掛器： 中國， 201020530991[P]. 20110330.

[10]曾鳴， 孫亮， 鐘朝廷， 等.海洋鋪管船用張緊器內(nèi)懸架的設計與分析[J].石油礦場機械， 2010， 39（6）： 3235.

ZENG Ming， SUN Liang， ZHONG Zhaoting， et al. Design and analysis in suspension of pipe tensioner machine used in offshore pipelay vessel[J]. Oil Field Equipment， 2010， 39（6）： 3235.

[11]朱江， 周文勝.中國近海油氣田區(qū)域開發(fā)戰(zhàn)略思考[J].中國海上油氣， 2009， 21（6）： 380382.

ZHU Jiang， ZHOU Wensheng. Pondering over regional development strategy for oil and gas fields offshore China[J]. China Offshore Oil & Gas， 2009， 21（6）： 380382.

[12]曾憲錦. 海上油氣田生產(chǎn)系統(tǒng)[M]. 北京： 石油工業(yè)出版社， 1993： 139156.

[13]李萬平. 計算流體力學[M]. 武漢： 華中科技大學出版社， 2004： 5458.

[14]劉玉林. 管道沿程水力損失的試驗研究[J]. 排灌機械， 1988， 6（5）： 2326.

LIU Yulin. Experimental study on hydraulic loss along the course of the pipeline[J]. Drainage & Irrigation Machinery， 1988， 6（5）： 2326.（編輯武曉英）endprint