陳世利，高春倩，郭世旭，黃新敬

CHEN Shili,GAO Chunqian,GUO Shixu,HUANG Xinjing

天津大學(xué) 精密測(cè)試技術(shù)及儀器國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津300072

State Key Laboratory of Precision Measuring Technology and Instruments,Tianjin University,Tianjin 300072,China

1 引言

海底管道作為海上油氣集輸?shù)闹饕侄我殉蔀楹Ｑ笥蜌赓Y源開發(fā)和運(yùn)輸?shù)纳€。由于海底表層地基不穩(wěn)定、介質(zhì)腐蝕、海流沖淘及海上作業(yè)等原因，海底油氣管道易發(fā)生油氣泄漏事故，造成巨大的經(jīng)濟(jì)損失與海洋環(huán)境污染[1-2]。為了消除管線運(yùn)行的潛在安全隱患，必須及時(shí)準(zhǔn)確地檢測(cè)和定位海底管道微小泄漏。

目前管道泄漏檢測(cè)方法可分為外檢測(cè)法和內(nèi)檢測(cè)法。廣泛應(yīng)用于陸地管道的負(fù)壓波法[3-4]、聲波相關(guān)法[5]、光纖泄漏法[6-8]等管道外檢測(cè)法僅能夠檢測(cè)大于總流量1%的原油泄漏和大于總流量5%的天然氣泄漏[9]，在海底管道的泄漏檢測(cè)中尚未見有效應(yīng)用，因?yàn)槿绱舜蟮挠蜌庑孤┰诤Ｃ嫔弦涯馨l(fā)現(xiàn)。海底油氣管道需要的是微小泄漏檢測(cè)技術(shù)。管道內(nèi)檢測(cè)法依靠?jī)?nèi)檢測(cè)器從管道內(nèi)靠近泄漏點(diǎn)監(jiān)聽泄漏信號(hào)，理論上能夠檢測(cè)到十分微小的油氣泄漏。傳統(tǒng)的管道內(nèi)檢測(cè)器（Pipeline Inspection Gauge，PIG）在前后壓力差的作用下前進(jìn)[10-12]，獲得管道內(nèi)壁的腐蝕、缺陷以及焊縫狀況等，通過壁面缺陷來判斷是否存在泄漏。但PIG 體積龐大，且與管壁緊密貼合，對(duì)管道變形十分敏感，存在高卡堵風(fēng)險(xiǎn)。而海底管道在洋流作用下易發(fā)生幾何變形，進(jìn)一步限制了其在海底管道的應(yīng)用。

為了解決海底管道微小泄漏檢測(cè)的迫切需求，本課題組于2007 年首次提出了一種球形內(nèi)檢測(cè)器[13]，能夠識(shí)別和精確定位0.1 L/min 的微小泄漏。球形載體在管道內(nèi)流體的推動(dòng)下滾動(dòng)前進(jìn)，沿途記錄管道內(nèi)泄漏信號(hào)和定位信息。由于球形載體外徑小于管道內(nèi)徑，且是滾動(dòng)前進(jìn)，所以不僅結(jié)構(gòu)簡(jiǎn)單，不易卡堵，而且噪聲小，檢測(cè)精度高[14-16]。然而，海底管道經(jīng)常存在高達(dá)數(shù)百米的立管段[17]，若球形內(nèi)檢測(cè)器不能順利通過立管段不僅無法實(shí)現(xiàn)微小泄漏檢測(cè)與定位，還會(huì)影響管道正常運(yùn)行，甚至誘發(fā)重大安全事故。因此對(duì)球形內(nèi)檢測(cè)器在海底管道立管段的通過能力的研究至關(guān)重要。

本文采用有限元仿真方法對(duì)球形內(nèi)檢測(cè)器在海底管道立管段的通過性進(jìn)行研究。為了驗(yàn)證仿真方法的正確性，搭建實(shí)驗(yàn)平臺(tái)管道，并設(shè)計(jì)不同模型球（球體密度、球體大小）。文中使用流體計(jì)算軟件Fluent 14.0，對(duì)實(shí)驗(yàn)平臺(tái)立管道內(nèi)模型球周圍流場(chǎng)進(jìn)行三維數(shù)值仿真并計(jì)算模型球順利通過立管段的臨界速度，通過與實(shí)驗(yàn)所得臨界速度進(jìn)行比較，驗(yàn)證仿真方法的正確性。采用該仿真方法，對(duì)實(shí)際海底管道立管段球形內(nèi)檢測(cè)器的通過性及影響因素進(jìn)行預(yù)測(cè)分析。

2 模型和數(shù)值分析

2.1 立管內(nèi)球體受力分析

以立管段內(nèi)球體為研究對(duì)象，如圖1 所示，理想條件下，球體在管道內(nèi)僅受到Ff,Ft和G的作用，F(xiàn)t為豎直向上的推力，由流體壓力以及流體速度能量作用于球體上產(chǎn)生，F(xiàn)f和G分別為球體所受浮力和重力是為已知量，豎直方向受力分析得：

圖1 豎直管道球體受力分析

顯然，當(dāng)Fh＞0 時(shí)，球體能順利通過豎直管道，使Fh=0 的Ft即為臨界推力，對(duì)應(yīng)流體速度為臨界流速。因此，求得Ft根據(jù)式（1），即可判斷球體是否能順利通過豎直管道，本文采用有限元仿真方法求解管內(nèi)流體對(duì)球體的推力Ft。

2.2 管內(nèi)流體數(shù)學(xué)模型

流體繞過球體在豎直管道內(nèi)的流動(dòng)為湍流，且流體為不可壓縮介質(zhì)，因此管道內(nèi)流體流動(dòng)遵循湍流運(yùn)輸方程，本文選用標(biāo)準(zhǔn)k-ε湍流模型。任何流體的流動(dòng)都滿足質(zhì)量守恒方程和動(dòng)量守恒方程，則笛卡爾坐標(biāo)系下管道內(nèi)流體流動(dòng)遵循如下控制方程[18]：

（1）質(zhì)量守恒方程

式中ρ表示管內(nèi)流體密度，t表示時(shí)間，u表示速度矢量，由于管內(nèi)流體流動(dòng)處于穩(wěn)態(tài)，則密度ρ不隨時(shí)間變化，公式（2）簡(jiǎn)化為：

（2）動(dòng)量守恒方程

式中，μ是動(dòng)力粘度，p是管內(nèi)流體微元體上的壓力，ui是u沿x、y、z方向的分速度，Si是動(dòng)量守恒方程的廣義源項(xiàng)，F(xiàn)i是微元體上的體力，由于管道內(nèi)流體是粘性為常數(shù)的不可壓縮流體，則si=0。

（3）不可壓縮流體的運(yùn)輸方程

其中，C1ε=1.44,C2ε=1.92,σk=1.0,σε=1.3,σk、σε分別是與湍動(dòng)能k和耗散率ε對(duì)應(yīng)的Prandtl 數(shù)，C1ε、C2ε為經(jīng)驗(yàn)常數(shù)。由以上控制方程經(jīng)流體計(jì)算軟件fluent14.0計(jì)算可得管道內(nèi)各微元體上的壓力p和速度u。

2.3 數(shù)值分析

由于直接搭建與實(shí)際管道尺寸相同的推力驗(yàn)證實(shí)驗(yàn)平臺(tái)十分困難，本文并未直接對(duì)球形內(nèi)檢測(cè)器在海底管道立管段的通過性進(jìn)行仿真計(jì)算和實(shí)驗(yàn)驗(yàn)證，而是將實(shí)際模型尺寸等比例縮小后進(jìn)行數(shù)值模擬，并針對(duì)小模型搭建實(shí)驗(yàn)平臺(tái)，設(shè)計(jì)和制作模型球，以驗(yàn)證仿真方法的正確性。實(shí)際球形內(nèi)檢測(cè)器樣機(jī)如圖2 所示，球殼密度為2 700 kg/m3的鋁殼，內(nèi)部搭載有加速度計(jì)、磁力計(jì)以及音波傳感器等，其平均密度為1 300 kg/m3。仿真計(jì)算和實(shí)驗(yàn)中對(duì)該極限密度（2 700 kg/m3）和平均密度（1 300 kg/m3）模型球在管道內(nèi)的運(yùn)動(dòng)進(jìn)行考察，根據(jù)實(shí)際選材難易情況，選擇密度為1 420 kg/m3聚甲醛和密度為2 700 kg/m3鋁分別制作模型球進(jìn)行研究。同時(shí)為了便于仿真結(jié)果和實(shí)驗(yàn)結(jié)果多組比較，設(shè)計(jì)多種尺寸實(shí)驗(yàn)?zāi)Ｐ颓颍睆絛分別為8 mm、16 mm、20 mm、24 mm、30 mm、32 mm。

圖2 實(shí)際球形內(nèi)檢測(cè)器

2.3.1 網(wǎng)格劃分

幾何模型如圖3 所示，管道內(nèi)徑D為40 mm，流體經(jīng)過模型球時(shí)，管內(nèi)流場(chǎng)經(jīng)歷了穩(wěn)定，劇烈變化到穩(wěn)定的過程，為了保證出入口處的流動(dòng)均是穩(wěn)定的，計(jì)算域上游L2 取1.5 m，下游L1 取4 m。采用ANSYS ICEM CFD 14.0 網(wǎng)格生成軟件對(duì)模型進(jìn)行網(wǎng)格劃分，生成網(wǎng)格如圖4。

圖3 幾何模型圖

圖4 六面體網(wǎng)格圖

整個(gè)計(jì)算域?yàn)榱骟w結(jié)構(gòu)化網(wǎng)格，相對(duì)于其他類型網(wǎng)格，使用更少的網(wǎng)格節(jié)點(diǎn)獲得更高精度結(jié)果。其中采用了3 個(gè)O-grid 對(duì)模型進(jìn)行細(xì)分，球面邊界外部用O-grid 進(jìn)行網(wǎng)格加密，采用等比例拉伸的結(jié)構(gòu)網(wǎng)格來加密球體附近管道網(wǎng)格，疏化遠(yuǎn)離球體管道網(wǎng)格，使網(wǎng)格能體現(xiàn)計(jì)算域內(nèi)流場(chǎng)的變化趨勢(shì)，網(wǎng)格數(shù)量約180 萬，經(jīng)檢驗(yàn)符合網(wǎng)格獨(dú)立性要求。

2.3.2 邊界條件和離散格式

實(shí)際應(yīng)用中流速是已知可調(diào)節(jié)量，因此入口邊界條件定義為速度入口，方向垂直于邊界，入口流速為0.2～1.5 m/s，出口邊界條件為壓力出口，其余邊界條件均定義為壁面邊界條件，默認(rèn)無滑移壁面條件，粗糙度常數(shù)為0.5，管內(nèi)流體為水，密度為999.4 kg/m3，動(dòng)力粘度為0.001 Pa·s，溫度以實(shí)驗(yàn)溫度為準(zhǔn)，為19 ℃。流體在管道內(nèi)流動(dòng)時(shí)，管壁區(qū)流場(chǎng)情況變化很大且存在粘性底層，標(biāo)準(zhǔn)k-ε湍流模型不能完全真實(shí)地描述管內(nèi)流體流動(dòng)，因此采用標(biāo)準(zhǔn)壁面函數(shù)改善近壁區(qū)域數(shù)值計(jì)算結(jié)果。流場(chǎng)計(jì)算采用Simple 算法進(jìn)行壓力修正，且用二階迎風(fēng)差分格式離散動(dòng)量方程。

2.3.3 仿真結(jié)果分析

采用前述數(shù)學(xué)模型、上述邊界條件和離散格式等，應(yīng)用Fluent 14.0軟件對(duì)管道內(nèi)流場(chǎng)進(jìn)行數(shù)值計(jì)算，圖5（a）和圖5（b）給出了流速相同時(shí)不同直徑球體所在管道流場(chǎng)的總壓力分布。

圖5 （a）直徑32 mm 球體總壓力云圖

圖5 （b）直徑16 mm 球體總壓力云圖

圖5（a）中球體直徑較大，球體上部?jī)蓚?cè)流場(chǎng)存在負(fù)壓區(qū)，上部附近流場(chǎng)為低壓區(qū)，這是由流體流經(jīng)球體后形成的尾跡發(fā)展的；圖5（b）中球體直徑較小，球體周圍流場(chǎng)的壓力均為正壓力，上部附近流場(chǎng)也存在低壓區(qū)，但面積較小。由圖5（a）和圖5（b）平面板的壓力值可以看出，直徑較大球體周圍流場(chǎng)壓力值遠(yuǎn)遠(yuǎn)高于直徑較小球體周圍流場(chǎng)的壓力值，由此可以初步斷定直徑大球體所受流體推力較大，相同密度球體通過同一管徑管道時(shí)所需流速較小。

3 仿真與實(shí)驗(yàn)對(duì)比

3.1 仿真結(jié)果計(jì)算

如圖5（a）和圖5（b），仿真計(jì)算出管內(nèi)流場(chǎng)各處壓力，應(yīng)用Fluent 計(jì)算流體對(duì)球體豎直方向推力。推力產(chǎn)生原因除壓力外，還包括粘性力。Fluent 計(jì)算球面沿豎直方向所受的合力大小等于球面各處壓力和粘性力與豎直向上矢量a點(diǎn)積加和的積分，計(jì)算公式如下：

其中：a·Fp為球面某點(diǎn)壓力沿豎直向上方向的分力，a·Fu為球面某點(diǎn)粘性力沿豎直向上方向的分力，a為沿豎直向上方向的單位矢量，F(xiàn)p為球面某點(diǎn)壓力矢量，F(xiàn)u為球面某點(diǎn)粘性力矢量。

通過上述方法計(jì)算得豎直管道內(nèi)流體對(duì)球體豎直方向的推力。本文采用有限元仿真方法對(duì)管道內(nèi)流體進(jìn)行仿真分析，由式（8）計(jì)算得Ft。為了得出不同模型球通過豎直管道的臨界速度，對(duì)仿真所得數(shù)據(jù)進(jìn)行處理分析。

圖6（a）為不同直徑模型球的推力-流速數(shù)據(jù)點(diǎn)曲線，觀察數(shù)據(jù)點(diǎn)曲線和二次多項(xiàng)式曲線極為相似，因此采用Polyfit[19]函數(shù)對(duì)推力和流速進(jìn)行二次多項(xiàng)式擬合，擬合曲線如圖6（b）所示，同時(shí)得到推力-流速擬合方程和數(shù)據(jù)點(diǎn)的相對(duì)誤差δ如表1 所示。

圖6 （a）不同直徑模型球推力-流速數(shù)據(jù)點(diǎn)曲線

圖6 （b）不同直徑模型球推力-流速擬合曲線

表1 推力-流速擬合方程

由表1 可知相對(duì)誤差在0.5%內(nèi)，可見擬合效果很好，推力與流速成二次函數(shù)關(guān)系。

已知模型球密度和直徑，結(jié)合公式（1）和擬合得到的推力-流速擬合方程，求得不同直徑聚甲醛和鋁質(zhì)模型球在管道內(nèi)能夠上升的臨界速度ul,如表2 所示。

表2 聚甲醛和鋁模型球的臨界速度

3.2 實(shí)驗(yàn)結(jié)果分析

圖7 為搭建的實(shí)驗(yàn)平臺(tái)流程圖，設(shè)計(jì)時(shí)主要考慮因素有：（1）為方便觀察模型球在管道內(nèi)的運(yùn)動(dòng)情況，管道選用透明有機(jī)玻璃管。但受有機(jī)玻璃加工工藝限制，彎管段選用鋼管，豎直和水平管道段為有機(jī)玻璃管。（2）考慮真實(shí)鋼材和鋼材彎頭的標(biāo)準(zhǔn)件尺寸，鋼管彎管段角度為90°和135°。（3）為方便實(shí)驗(yàn)時(shí)能多次快速地投球和收球，設(shè)計(jì)了收發(fā)球裝置。（4）在流體入口處安裝有針型閥和流量計(jì)，用于調(diào)節(jié)和顯示管道內(nèi)流體流速。（5）為避免污染實(shí)驗(yàn)場(chǎng)地，實(shí)驗(yàn)裝置整體組成一個(gè)循環(huán)系統(tǒng)，流體在管道內(nèi)循環(huán)流動(dòng)。

圖7 實(shí)驗(yàn)平臺(tái)流程圖

管道內(nèi)流體為水，動(dòng)力粘度為0.001 Pa·s，溫度為19 ℃。實(shí)驗(yàn)時(shí)水池中的水經(jīng)水泵流入實(shí)驗(yàn)管道內(nèi)，當(dāng)水充滿管道后，將不同直徑的聚甲醛和鋁質(zhì)模型球分別投入實(shí)驗(yàn)管道，調(diào)節(jié)針型閥改變管道內(nèi)水的流速使球體在豎直管道內(nèi)懸浮，通過流量計(jì)讀取其臨界速度ul',最后增大流速對(duì)模型球進(jìn)行回收。表3 為實(shí)驗(yàn)讀取的臨界速度ul'及實(shí)驗(yàn)與仿真結(jié)果的誤差e。

實(shí)驗(yàn)與仿真存在一定誤差，誤差產(chǎn)生主要有以下原因：

（1）針型閥調(diào)節(jié)流速時(shí)流量計(jì)讀數(shù)反應(yīng)慢，與管內(nèi)流體流速存在一定滯后，流量計(jì)讀數(shù)有一定偏差。

（2）實(shí)驗(yàn)條件與仿真條件存在一定的差異，如調(diào)節(jié)流速時(shí)存在壓力的突變等。

表3 實(shí)驗(yàn)結(jié)果及誤差

（3）流體中夾雜少量氣泡，影響實(shí)驗(yàn)結(jié)果。

由表3 可知，誤差范圍為－0.05～0.05 m/s，此誤差在工程應(yīng)用中是可接受的，仿真結(jié)果和實(shí)驗(yàn)結(jié)果基本吻合。這說明此仿真方法是正確的，仿真結(jié)果真實(shí)可靠，可以用于預(yù)測(cè)球形內(nèi)檢測(cè)器海底管道立管段的通過能力。

4 通過性影響因素分析

應(yīng)用上述仿真方法分別對(duì)球形內(nèi)檢測(cè)器在8～28 inch立管段時(shí)管內(nèi)流場(chǎng)進(jìn)行數(shù)值模擬，實(shí)際工況下管內(nèi)流體流速為0.8～1.2 m/s，因此選取仿真流速為0.8～1.2 m/s，流體為石油。結(jié)合公式（6）計(jì)算流體對(duì)球形內(nèi)檢測(cè)器的推力，再結(jié)合公式（1）得出球形內(nèi)檢測(cè)器在海底管道立管段受到的合力，應(yīng)用差值方法[20-21]對(duì)合力做曲面擬合，擬合曲面如圖8 所示。

圖8 實(shí)際工況合力擬合曲面

由圖8可以看出，當(dāng)管徑不變時(shí)，隨著流速的增加球形內(nèi)檢測(cè)器的合力增大，但當(dāng)管道直徑較大時(shí)推力隨流速增加不明顯，且合力均小于零，球形內(nèi)檢測(cè)器不能通過海底管道立管段；當(dāng)管道直徑為10 in 和8 in 時(shí)，即球管徑比達(dá)到73%時(shí)，其合力明顯增大，且合力隨流速增加而增加顯著，此時(shí)合力均大于零，球形內(nèi)檢測(cè)器均能通過海底管道立管段?？傻媒Y(jié)論：球形內(nèi)檢測(cè)器直徑不變時(shí)，管道直徑越小，即球管徑比越大，球形內(nèi)檢測(cè)器所受推力和合力越大。球管徑比是一個(gè)重要的參數(shù)，也是影響球體在管道立管段通過性的主要因素，當(dāng)球管徑比達(dá)到73%時(shí)，球形內(nèi)檢測(cè)器在正常工作流速0.8～1.2 m/s內(nèi)，均能順利通過海底管道立管段。

5 結(jié)束語

本文對(duì)采用有限元仿真方法對(duì)球形內(nèi)檢測(cè)器安全通過性及其影響因素進(jìn)行研究，主要得到了以下結(jié)論。

（1）通過Fluent 仿真立管道內(nèi)球體周圍三維流場(chǎng)，比較仿真結(jié)果和實(shí)驗(yàn)結(jié)果，驗(yàn)證了仿真方法是正確的，其結(jié)果真實(shí)可信，此仿真方法可用于預(yù)測(cè)對(duì)球形內(nèi)檢測(cè)器在海底管道立管段的通過性。

（2）在管徑和球體直徑保持不變時(shí)，豎直管道內(nèi)球體所受推力和流速呈二次函數(shù)關(guān)系。

（3）球管徑比對(duì)球形內(nèi)檢測(cè)器在海底管道立管內(nèi)的通過性來說是一個(gè)很重要的參數(shù)，球管徑比大于等于73%時(shí)，在其正常工況流速（0.8～1.2 m/s）下即可通過。

[1] 覃覓.石油泄漏災(zāi)難[J].資源與人居環(huán)境，2006（11S）：64-67.

[2] 曉武.石油泄漏，海洋面臨生態(tài)之痛[J].節(jié)能與環(huán)保，2011（9）：36-39.

[3] 邵煜，葛傳虎，葉昊，等.基于負(fù)壓波的管道泄漏檢測(cè)與定位系統(tǒng)評(píng)價(jià)[J].油氣儲(chǔ)運(yùn)，2008，27（4）：5-9.

[4] Tian Chunhua，Yan Junchi，Huang Jin，et al.Negative pressure wave based pipeline leak detection：challenges and algorithms[C]//Proceedings of IEEE International Conference on Service Operations and Logistics and Informatics，2012.

[5] 華科，葉昊，王桂增，等.基于聲波的輸氣管道泄漏檢測(cè)與定位技術(shù)[J].華中科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，37（8）：181-183.

[6] Tu Y，Chen H.Design of oil pipeline leak detection and communication system based on optical fiber technology[C]//Proceedings of SPIE 1999，1999，3737：584-592.

[7] 周琰，靳世久，張昀超，等.分布式光纖管道泄漏檢測(cè)和定位技術(shù)[J].石油學(xué)報(bào)，2006，27（2）：121-124.

[8] Paranjape R，Liu N，Rumple C，et al.A distributed fiber optic system for oil pipeline leakage detection[C]//Proc of SPIE 2003，F(xiàn)eb 1，2003，4833.

[9] 袁朝慶，龐鑫峰，劉燕.管道泄漏檢測(cè)技術(shù)現(xiàn)狀及展望[J].大慶石油學(xué)院學(xué)報(bào)，2006（2）：76-79.

[10] 宋志東，靳世久，李一博，等.管道內(nèi)檢測(cè)器里程輪信號(hào)優(yōu)選算法的設(shè)計(jì)與實(shí)現(xiàn)[J].管道技術(shù)與設(shè)備，2006（3）：12-13.

[11] 劉保余，綦耀光.輸氣管道內(nèi)檢測(cè)器設(shè)計(jì)及理論研究[D].北京：中國(guó)石油大學(xué)，2010.

[12] 張立軍，綦耀光，劉保余，等.輸氣管道內(nèi)檢測(cè)器模型試驗(yàn)系統(tǒng)的設(shè)計(jì)[J].管道技術(shù)與設(shè)備，2010（2）：26-27.

[13] 郭世旭.檢測(cè)輸油管道泄漏的球形內(nèi)檢測(cè)器：中國(guó)，201110402319.3[P].2012-05-09.

[14] Chen Shili，Guo Tingting.A length measuring method based on accelerometer[C]//Proceedings of the 1st International Conference on Instrumentation，Measurement，Computer，Communication and Control，2003.

[15] 劉淵，陳世利，郭世旭，等.輸油管道微小泄漏球形檢測(cè)器的設(shè)計(jì)與實(shí)現(xiàn)[J].傳感器與微系統(tǒng)，2012，31（6）：130-133.

[16] 趙偉.油氣管道內(nèi)磁場(chǎng)特性分析及應(yīng)用研究[D].天津：天津大學(xué)，2012.

[17] 黃鑫，余建興.海底管線設(shè)計(jì)分析[D].天津：天津大學(xué)，2004.

[18] 王福軍.計(jì)算流體動(dòng)力學(xué)分析[M].北京：清華大學(xué)出版社，2004.

[19] 饒秀勤，應(yīng)義斌.自適應(yīng)濾波器在水果聲學(xué)特性研究中的應(yīng)用[J].浙江大學(xué)學(xué)報(bào)：農(nóng)業(yè)與生命科學(xué)版，2003（2）：28-31.

[20] 彭芳瑜，周云飛，周濟(jì).基于插值與逼近的復(fù)雜曲面擬合[J].工程圖學(xué)學(xué)報(bào)，2002（4）：88-96.

[21] 蘭鳳崇，陳吉清，鄭世紅，等.散亂數(shù)據(jù)的曲面擬合及其在汽車車身外表面造型設(shè)計(jì)中的應(yīng)用[J].吉林工業(yè)大學(xué)學(xué)報(bào)，1996（3）：10-15.