陳 靖，朱克強(qiáng)

(寧波大學(xué) 海運(yùn)學(xué)院， 浙江 寧波 315211)

0 引 言

近海平臺(tái)和海底管系的各結(jié)構(gòu)部件經(jīng)常暴露在海水中，分析和計(jì)算波浪對海洋結(jié)構(gòu)物的作用力是進(jìn)行此類海洋結(jié)構(gòu)物設(shè)計(jì)建造的基礎(chǔ)。當(dāng)管件或者結(jié)構(gòu)部件的直徑與該處波浪波長的比率小于0.2時(shí)，我們稱此類海洋結(jié)構(gòu)物為“小構(gòu)件”，小構(gòu)件在海水中受力情況通常利用莫里森方程來計(jì)算[1]，莫里森方程的微分形式如式(1)。影響莫里森方程計(jì)算結(jié)果的一個(gè)主要因素是方程中慣性力系數(shù)和阻力系數(shù)的取值，針對小構(gòu)件圓柱體而言，對于慣性力系數(shù)和阻力系數(shù)的取值，通常做法是根據(jù)表1波浪力載荷計(jì)算方法指南來取值。這種取值方式雖然將取值過程簡單化，但是該表系數(shù)取值的前提條件較多，取得系數(shù)值浮動(dòng)性較大，很容易導(dǎo)致最后計(jì)算誤差較大。

目前，針對慣性力系數(shù)與阻力系數(shù)取值的相關(guān)研究很多，比較典型的有Brown和Lawler[2](見式(2))，以及通過各種數(shù)學(xué)和實(shí)驗(yàn)方法對此公式具體系數(shù)進(jìn)一步改進(jìn)[3-4]。Pilyguin等[5]研究了雷諾數(shù)在4×105～1×106范圍內(nèi)，流體不同流速對慣性力系數(shù)取值的影響。Koji Otsuka等[6]在分析垂直圓柱體在低邱卡數(shù)下慣性力受力情況時(shí)，得到慣性力系數(shù)與邱卡數(shù)的大致函數(shù)關(guān)系如式(3)。大部分的研究主要通過歸納試驗(yàn)結(jié)果并對結(jié)果進(jìn)行數(shù)學(xué)推算，由于研究的針對性，很少將邱卡數(shù)和雷諾數(shù)同時(shí)列入考慮。本文考慮利用數(shù)學(xué)軟件對原始試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)進(jìn)行曲線擬合，得到圓柱體慣性力系數(shù)和阻力系數(shù)隨雷諾數(shù)和邱卡數(shù)變化曲線的函數(shù)。

dF=CmρdVu·n+CdρdAu0u0， (1)

1 離散點(diǎn)組采集和曲線擬合

1.1 離散點(diǎn)組采集

所謂擬合是指平面上離散點(diǎn)組所表示的坐標(biāo)之間的函數(shù)關(guān)系的一種數(shù)據(jù)處理方法。從試驗(yàn)與經(jīng)驗(yàn)中得到若干離散函數(shù)值{f1,f2,…,fn}，以這些給定的數(shù)據(jù)點(diǎn)為基礎(chǔ)，通過調(diào)整該函數(shù)中若干待定系數(shù)，確定滿足特定要求的函數(shù)曲線，并使該函數(shù)與已知點(diǎn)集的差別盡量小，稱之為曲線擬合[7]。

本文根據(jù)試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)[8]，以雷諾數(shù)Re為笛卡兒平面上X軸，以慣性力系數(shù)(Cm)為Y軸，由于篇幅原因此處以邱卡數(shù)等于40時(shí)，圓柱體的慣性力系數(shù)隨雷諾數(shù)變化為例，初步采集X和Y的離散點(diǎn)數(shù)據(jù)如下：Y=[1.02 1.05 1.1 1.16 1.25 1.37 1.42 1.45 1.5 1.55 1.6 1.69 1.77 1.81 1.82 1.81 1.78]；X=[0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.25 1.5 2 2.5 2.75 3 4]。初步采集的離散點(diǎn)組可能擁有較大誤差，因此必須進(jìn)行修正。將初步采集的數(shù)據(jù)輸入Matlab，生成曲線圖象比對試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)圖象，對個(gè)別離散點(diǎn)值進(jìn)行盡量接近試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)圖象的數(shù)據(jù)修正，得到的結(jié)果如下：Y=[1.02 1.05 1.1 1.16 1.26 1.33 1.41 1.45 1.5 1.55 1.6 1.69 1.77 1.81 1.82 1.81 1.78]；X=[0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.25 1.5 2 2.5 2.75 3 4]。這里需要說明的是，進(jìn)行最終擬合的離散點(diǎn)組并不需要完全符合試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)有2個(gè)主要原因：1)由于試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)也并非精確值，修正過的離散點(diǎn)數(shù)據(jù)誤差很小，因此足以用于實(shí)踐計(jì)算；2)曲線擬合的最終結(jié)果函數(shù)也并不會(huì)準(zhǔn)確經(jīng)過所有離散點(diǎn)組，因此過于精確的數(shù)據(jù)對于最終的擬合結(jié)果意義不大。

1.2 曲線擬合的基本原理

在現(xiàn)實(shí)的物理研究過程中，常常需要得到多個(gè)測量數(shù)據(jù)變量的函數(shù)關(guān)系，此時(shí)通常使用一元曲線擬合方法，將實(shí)驗(yàn)或經(jīng)驗(yàn)數(shù)據(jù)結(jié)合數(shù)學(xué)方法得到多個(gè)物理量之間的近似函數(shù)表達(dá)式。多數(shù)情況下，選擇曲線使得數(shù)據(jù)點(diǎn)的平方誤差和最小。這種選擇就是最小二乘曲線擬合，即曲線擬合最常用的方法，其基本思路如式(4)，其中ki(x)是事先選定的一組函數(shù)；ai為待定系數(shù)(i=0,1,2…m,m

(4)

1.3 運(yùn)用Matlab進(jìn)行曲線擬合

1.3.1 曲線擬合工具Cftool簡介

Cftool是Matlab7.0內(nèi)置的曲線擬合工具箱，擁有可視化的圖形界面和強(qiáng)大的圖形擬合功能，包括：1)可視化地展開一個(gè)或多個(gè)數(shù)據(jù)集，并可用散點(diǎn)圖來表示；2)用殘差和置信區(qū)間可視化地估計(jì)擬合結(jié)果的好壞；3)通過其他界面可以實(shí)現(xiàn)許多附加功能：輸出、察看和平滑數(shù)據(jù)，擬合數(shù)據(jù)、比較擬合曲線和數(shù)據(jù)集，從擬合曲線中排除特殊的數(shù)據(jù)點(diǎn)，選定區(qū)間后可以顯示擬合曲線和數(shù)據(jù)集，還可以對數(shù)據(jù)進(jìn)行內(nèi)插法、外推法、微分或積分?jǐn)M合[9]。Cftool提供的擬合類型有：指數(shù)逼近、傅立葉逼近、高斯逼近、插值逼近、多項(xiàng)式逼近、正弦曲線逼近、冪逼近等，適合于各種復(fù)雜模型的曲線擬合。

1.3.2 曲線擬合具體操作

由于篇幅問題，無法對圓柱體的慣性力系數(shù)和阻力系數(shù)隨雷諾數(shù)和邱卡數(shù)變化曲線的擬合過程在文中全部展現(xiàn)，現(xiàn)仍以邱卡數(shù)等于40時(shí)，圓柱體的慣性力系數(shù)隨雷諾數(shù)變化為例，在Matlab中可以直接通過Cftool命令調(diào)用曲線擬合工具箱，對其進(jìn)行曲線擬合。

在Matlab主窗口輸入修正后離散點(diǎn)數(shù)據(jù)，運(yùn)用Cftool(X,Y)命令設(shè)置擬合數(shù)據(jù)坐標(biāo)系，并打開擬合工具箱，單擊“Fitting”進(jìn)行數(shù)據(jù)擬合設(shè)置。分別設(shè)置擬合多項(xiàng)式次數(shù)為3次、4次、5次和6次，具體擬合效果如圖1所示。由圖1可以得到5次多項(xiàng)式擬合時(shí)，曲線起伏較小且貼近離散點(diǎn)組，同時(shí)也符合試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)圖象，擬合效果最好。此處需要說明的是：由圖1可得，此例中當(dāng)雷諾數(shù)大于4×105時(shí)，慣性力系數(shù)發(fā)生的變化可以忽略不計(jì)，因此不加入Cftool擬合中。以此類推，重復(fù)數(shù)據(jù)采集和曲線擬合的操作，可以得到圓柱體慣性力系數(shù)和阻力系數(shù)隨雷諾數(shù)和邱卡數(shù)變化曲線的所有擬合結(jié)果。

圖1 離散點(diǎn)組不同多項(xiàng)式次數(shù)擬合結(jié)果對比Fig.1 The discrete point group fitting comparison result with different degrees of polynomial

1.3.3 擬合效果評價(jià)

對于擬合的效果評價(jià)，通常使用2種指標(biāo)：負(fù)相關(guān)系數(shù)與均方根誤差。前文例子中Matlab給出擬合評價(jià)指標(biāo)如下：SSE(誤差平方和):0.003 922；R-square(復(fù)相關(guān)系數(shù)或復(fù)測定系數(shù)):0.997；Adjusted R-square(調(diào)整自由度復(fù)相關(guān)系數(shù)):0.995 6；RMSE(均方根誤差):0.018 88。可以看出，SSE和RMSE指標(biāo)比較小，且R-square接近于1 ，說明擬合效果較好，具有使用價(jià)值。

1.4 擬合成果的執(zhí)行程序編寫

得到曲線擬合的結(jié)果之后，由于得到的擬合函數(shù)和相關(guān)數(shù)據(jù)較為繁多，需要一個(gè)簡單的可視化執(zhí)行程序，以便在進(jìn)行相關(guān)計(jì)算時(shí)，能夠準(zhǔn)確快速地得到慣性力系數(shù)和阻力系數(shù)。本文運(yùn)用Visual Basic 6.0將擬合結(jié)果編寫為一個(gè)exe執(zhí)行程序，具體編寫思路如下：首先建立若干控件，包括：用于輸入雷諾數(shù)的文本框(編寫相應(yīng)代碼，確保用戶輸入數(shù)據(jù)的正確性，減少程序運(yùn)行出錯(cuò)幾率)；用于選擇邱卡數(shù)值得下拉列表框；用于開啟執(zhí)行程序的命令按鈕；用于顯示擬合得到的慣性力和阻力系數(shù)結(jié)果以及其他說明的標(biāo)簽若干。執(zhí)行程序的核心代碼是擬合函數(shù)編寫與調(diào)用，仍然以前文擬合結(jié)果為例，定義擬合結(jié)果函數(shù)f40代碼如下：

Private Sub f40()

x=Text1; If x>4 Then Label1=1.77; ElseIf x>=0.15 And x<=4 Then Const p1=-0.004975

……

Const p6=0.8759; Label1=p1*x^5+p2*x^4+p3*x^3+p4*x^2+p5*x+p6

Else Label1= "無相應(yīng)數(shù)據(jù)";End If;If x>4 Then Label2=0.6

ElseIf x>=0.19 And x<=4 Then Const q1=0.01435

……

Const q8=1.583; Label2=q1*x^7+q2*x^6+q3*x^5+q4*x^4+q5*x^3+q6*x^2+q7*x+q8; Else Label2="無相應(yīng)數(shù)據(jù)"; End If; End Sub

選擇邱卡數(shù)后，調(diào)用對應(yīng)擬合函數(shù)的代碼如下：

Private Sub Command1_Click()

If Combo1.Text="6" Then Call f6; ElseIf Combo1.Text="8" Then Call f8

……

ElseIf Combo1.Text="100" Then Call f100; End If; End Sub。

2 算 例

通過一個(gè)算例來驗(yàn)證曲線擬合結(jié)果的實(shí)際意義。假設(shè)線性波理論適用，一個(gè)固定導(dǎo)管架平臺(tái)在30 m水深海域，該處波浪周期為11 s，此結(jié)構(gòu)物主要樁腿是垂直于海水的圓柱鋼管，直徑1 m。試計(jì)算單個(gè)樁腿在1個(gè)波浪周期內(nèi)的阻力、慣性力和合力的變化。

傳統(tǒng)算法計(jì)算結(jié)果如表2所示，其中慣性力系數(shù)和阻力系數(shù)根據(jù)表1來取值。慣性力系數(shù)和阻力系數(shù)按照曲線擬合結(jié)果根據(jù)具體邱卡數(shù)和雷諾數(shù)得到的計(jì)算結(jié)果如表3所示，其中邱卡數(shù)取近似值15。

表2 系數(shù)查表計(jì)算圓柱體波浪周期內(nèi)阻力、慣性力和合力的變化

表3 系數(shù)擬合取值計(jì)算圓柱體波浪周期內(nèi)阻力、慣性力和合力的變化

比對計(jì)算結(jié)果，由于2個(gè)重要系數(shù)取值的差異，最終計(jì)算結(jié)果也產(chǎn)生一定差異。圓柱體受到的最大慣性力和最大合力，擬合取值算法結(jié)果比查表取值算法減少了2.78%，誤差小于5%在可接受范圍內(nèi)。然而圓柱體所受最大阻力，擬合取值算法結(jié)果比查表取值算法增加25.81%，誤差大于5%。根據(jù)擬合結(jié)果得到慣性力系數(shù)和阻力系數(shù)相比查表取值更接近現(xiàn)實(shí)，得出的計(jì)算結(jié)果更具實(shí)踐指導(dǎo)意義。

3 結(jié) 語

莫里森方程在海洋工程領(lǐng)域應(yīng)用廣泛，方程中慣性力系數(shù)和阻力系數(shù)的取值直接影響最終工程計(jì)算結(jié)果。因此慣性力系數(shù)和阻力系數(shù)越精確，理論模型的計(jì)算結(jié)果就越具有實(shí)踐意義。本文利用Matlab軟件的Cftool擬合工具箱，對圓柱體慣性力系數(shù)和阻力系數(shù)隨雷諾數(shù)和邱卡數(shù)變化的試驗(yàn)與經(jīng)驗(yàn)數(shù)據(jù)進(jìn)行擬合，再編寫成相應(yīng)的執(zhí)行程序(具體編寫代碼略)，能夠方便快速地得到較為準(zhǔn)確的慣性力系數(shù)和阻力系數(shù)。通過一個(gè)簡單算例證明，擬合結(jié)果得到慣性力系數(shù)和阻力系數(shù)相比查表取值更接近實(shí)際情況，因此擬合有較強(qiáng)的實(shí)際意義。本文曲線擬合使用的是多項(xiàng)式擬合的方式，同時(shí)由于慣性力系數(shù)和阻力系數(shù)準(zhǔn)確取值需要考慮的因素很多，擬合結(jié)果取值的精確性還有待提高。

[1] RANDALL R E.Elements of ocean engineering[M].上海交通大學(xué)出版社,2002.

[2] BROWN P P,LAWLER D F.Sphere drag and setting velocity revisited[J], Journal of Environmental Engineering-ASCE,2003(3)：222-231.

[3] CLIF R,GRACE J R,WEBER M E.Bubbles,drops,and particles[M].Academic Press,New York,1978.

[4] TURTON R,LEVENSPIEL O.A short note on the drag correlation for spheres[J].Powder Technology,1986,47(1):83-86.

[5] PILYUGIN N N,KHLEBNIKOV V S.Investigation of aer-odynamic drag of two bodies in trans- and supersonic flows[J].Journal of Applied Mechanics and Technical Physics,2003,44(2):187-192.

[6] KOJI O,YOSHIHO I.Estimation of inertia forces on a horizontal circular cylinder in regular and irregular waves at low Keulegan-Carpenter numbers[J].Applied Ocean Research,1996,18:145-156.

[7] 樊曉紅.樣條函數(shù)與基于樣條函數(shù)的數(shù)據(jù)擬合方法[J].科技信息,2011(9)：507-508.

[8] SARPKAYA,ISAACSON.Mechanics of wave forces on offshore structures[J]，New York:Van Nostrand Reinhold Co.,1981.

[9] 蘇金明,張蓮花,劉波,等.Matlab工具箱應(yīng)用[M].北京:電子工業(yè)出版社,2004:489-512.

[10] 倫冠德.MATLAB曲線擬合工具箱在試驗(yàn)數(shù)據(jù)處理上的應(yīng)用[J].拖拉機(jī)與農(nóng)用運(yùn)輸車,2006(8)：90-91.