三體船興波阻力計(jì)算的自由面網(wǎng)格快速生成

王 中,盧曉平,王 瑋

(1.海軍工程大學(xué) 船舶與動(dòng)力學(xué)院,湖北 武漢 430033;2.海軍裝備研究院 艦船所,北京 100073)

采用三維Rankine源方法計(jì)算船舶興波阻力時(shí),自由面網(wǎng)格的生成是一個(gè)必要步驟,也是一個(gè)關(guān)鍵步驟[1].近年來(lái)三體船成為國(guó)內(nèi)外研究的熱點(diǎn)[2-4].由于側(cè)體的存在,自由面網(wǎng)格的劃分比單體船更加復(fù)雜、繁瑣.Dawson型方法自由面采用流線網(wǎng)格,三體船側(cè)體附近的流線需要反復(fù)迭代計(jì)算、調(diào)整才能確定,費(fèi)時(shí)費(fèi)力;另外,更為合理的自由面條件不能完整地用流線的導(dǎo)數(shù)表示,如改進(jìn)的Dawson法自由面條件,包括Rankine源方法各種階次非線性自由面條件以及全非線性自由面條件.在這一類(lèi)興波阻力計(jì)算方法中,舍棄流線網(wǎng)格,采用自由面上的“貼水線”網(wǎng)格勢(shì)在必行[5-6].但對(duì)三體船來(lái)說(shuō),快速生成自由面上的“貼水線”網(wǎng)格也并非易事.側(cè)體以及方尾的存在使得網(wǎng)格生成變得相當(dāng)復(fù)雜,采用商用軟件(如Gambit、ICEM等)來(lái)劃分此類(lèi)網(wǎng)格固然可以解決問(wèn)題,但過(guò)程繁瑣,此類(lèi)軟件生成的網(wǎng)格文件必須經(jīng)過(guò)進(jìn)一步的解析轉(zhuǎn)換才能供自己編制的興波阻力計(jì)算程序使用.因此,針對(duì)三體船構(gòu)型特點(diǎn)開(kāi)發(fā)適用于三體船的自由面網(wǎng)格快速生成系統(tǒng),成為提高三體船興波問(wèn)題計(jì)算效率的重要技術(shù)環(huán)節(jié).

1 網(wǎng)格生成方法概述

實(shí)用的自由面網(wǎng)格快速生成系統(tǒng)應(yīng)符合輸入簡(jiǎn)單、修改靈活、存儲(chǔ)方便、過(guò)程可視等要求,自由面網(wǎng)格快速生成系統(tǒng)界面如圖 1所示,具備上述特點(diǎn).通過(guò)該系統(tǒng),可將整個(gè)輸入?yún)?shù)及網(wǎng)格劃分結(jié)果保存于數(shù)據(jù)文件中,需要修改時(shí)直接打開(kāi)該文件,改變一些參數(shù)即可快速重新生成需要的網(wǎng)格;同時(shí),可將生成的網(wǎng)格數(shù)據(jù)導(dǎo)出成文本文件,直接供其他程序使用,導(dǎo)出的網(wǎng)格文件分塊存儲(chǔ)了主區(qū)域、主體方尾后側(cè)區(qū)域、側(cè)體方尾后側(cè)區(qū)域 3部分網(wǎng)格頂點(diǎn)坐標(biāo)信息.

在網(wǎng)格劃分前首先需要確定必要的信息.按如下方式定義系統(tǒng)的輸入?yún)?shù),如圖 1、2所示:主體靜水面水線型值點(diǎn)、側(cè)體靜水面水線型值點(diǎn)(側(cè)體局部坐標(biāo)系)及側(cè)體橫向偏距Dy和縱向偏距Dx、上游長(zhǎng)度P1P12、下游長(zhǎng)度P13P2、橫向?qū)挾?P1P11、傾角∠P11P1P8、縱向單位船長(zhǎng)網(wǎng)格數(shù)、橫向總網(wǎng)格數(shù)、主體與側(cè)體間橫向網(wǎng)格數(shù)、橫向網(wǎng)格密度控制因子、主體方尾橫向網(wǎng)格數(shù)、側(cè)體方尾橫向網(wǎng)格數(shù),其中上游長(zhǎng)度、下游長(zhǎng)度、橫向長(zhǎng)度以主體長(zhǎng)度 L為單位.

圖1 軟件系統(tǒng)主界面Fig.1 Window of the software

圖2 自由面網(wǎng)格參數(shù)定義Fig.2 Free-surfacemesh parameter definition

如圖 2所示,由輸入的網(wǎng)格控制參數(shù)可以確定邊界點(diǎn):P1、P2、P3、P4、P5、P6、P7、P8、P9、P10.這些邊界點(diǎn)與輸入的船體水線型值點(diǎn)可生成供插值的NURBS縱向邊界線:l1、l2、l3、l4、l5、l6、l7、l8、l9,參見(jiàn)圖2.這些邊界線是由一條或多條NURBS曲線組成的,不論是 2點(diǎn)組成的線段還是多點(diǎn)組成的曲線,均統(tǒng)一寫(xiě)成NURBS曲線形式,邊界線可以將自由面劃分成 S1～S5共 5個(gè)子區(qū)域,然后根據(jù)定義的縱向和橫向網(wǎng)格數(shù)在各個(gè)子區(qū)域分別劃分網(wǎng)格.NURBS曲線表示三體船自由面網(wǎng)格(子)區(qū)域的邊界線是本文貼體坐標(biāo)代數(shù)法生成自由面網(wǎng)格的特點(diǎn)之一,這種方法便于自由面網(wǎng)格邊界線分割,從而較好地控制網(wǎng)格線密度.

一旦按上述方法確定了自由面網(wǎng)格的邊界線,即可以分區(qū)采用貼體坐標(biāo)代數(shù)法生成自由面網(wǎng)格.自由面橫向網(wǎng)格和縱向網(wǎng)格密度的控制對(duì)興波阻力Rankine源方法計(jì)算結(jié)果的影響不可忽視,在貼體坐標(biāo)代數(shù)法生成自由面網(wǎng)格過(guò)程中可采用各種方法控制網(wǎng)格密度,本文縱向網(wǎng)格(即橫向網(wǎng)格線)取為等間距的,重點(diǎn)考慮了橫向網(wǎng)格(即縱向網(wǎng)格線)密度的控制.經(jīng)分析和試算,確定在 l4邊界以內(nèi)采用較為密集的等間距橫向網(wǎng)格,而在 l4邊界以外的 S3區(qū)域采用等比數(shù)列方法控制橫向網(wǎng)格密度,假設(shè)對(duì)一段NURBS曲線,要按照等比因子k分成N份,即將參數(shù)域u∈[0,1]按k分成N份,然后根據(jù)得到的u0、u1、u2、…、uN計(jì)算NURBS曲線上的點(diǎn)P(ui),根據(jù)得到的這些點(diǎn)即可得到網(wǎng)格的邊界線段,如圖 3所示.

圖3 NURBS曲線分段示意圖Fig.3 NURBS curve subsection

設(shè)第i個(gè)分段參數(shù)變化量:

則有

由式(1)可以計(jì)算得到

以上橫向網(wǎng)格(即縱向網(wǎng)格線)密度的控制是本文貼體坐標(biāo)代數(shù)法生成自由面網(wǎng)格的特點(diǎn)之二.歸納以上自由面網(wǎng)格生成方法,可以給出本系統(tǒng)劃分網(wǎng)格的流程如圖4所示.

圖 4 自由面網(wǎng)格快速生成流程圖Fig.4 Flow chart of free-surface mesh generation

以上探討的是方尾三體船型自由面網(wǎng)格的劃分,而其他船型,如不含方尾的三體、雙體、單體船型,自由面網(wǎng)格劃分均屬于它的子集,更容易生成.

2 網(wǎng)格生成示例

圖 5所示為采用該系統(tǒng)生成的自由面網(wǎng)格,其中圖5(a)是某單體方尾船型自由面網(wǎng)格劃分軟件屏幕截圖,自由面計(jì)算區(qū)域參數(shù)設(shè)置與圖 1相同,縱向單位船長(zhǎng)網(wǎng)格數(shù)為 20,橫向總網(wǎng)格數(shù)為 20,橫向網(wǎng)格密度控制因子為 1.07,方尾橫向劃分 2個(gè)網(wǎng)格,總網(wǎng)格數(shù)為1 260.

圖5(b)給出的是某5 000 t級(jí)三體船自由面網(wǎng)格劃分截圖,自由面計(jì)算區(qū)域同圖1參數(shù),縱向單位船長(zhǎng)網(wǎng)格數(shù)為 30,橫向總網(wǎng)格數(shù)為 30,主體與側(cè)體間橫向網(wǎng)格數(shù)為 8,傾角為 45°,橫向網(wǎng)格密度控制因子為 1.1,半個(gè)主體方尾橫向網(wǎng)格數(shù)為 3,側(cè)體方尾橫向網(wǎng)格數(shù)為 2,總網(wǎng)格數(shù)為 2 935.由圖 5可以看出,所開(kāi)發(fā)的自由面網(wǎng)格快速生成系統(tǒng)能夠快速生成較復(fù)雜的自由面網(wǎng)格,簡(jiǎn)化了網(wǎng)格劃分的過(guò)程;而且通過(guò)合理控制橫向網(wǎng)格密度,節(jié)省了自由面單元數(shù),提高了自由面網(wǎng)格生成效率.

圖 5 自由面網(wǎng)格劃分示例截圖Fig.5 Free-surfacemesh examp les

3 數(shù)值實(shí)驗(yàn)

采用以上自由面網(wǎng)格生成方法基于Wigley數(shù)學(xué)三體船[8]實(shí)施興波阻力數(shù)值計(jì)算,并對(duì)計(jì)算結(jié)果進(jìn)行分析.Wigley數(shù)學(xué)三體船主體和側(cè)體船型函數(shù)表達(dá)式在各自坐標(biāo)系下均為

式中:L/B=10,B/T=1.6,側(cè)體長(zhǎng)度為中體長(zhǎng)度的一半.根據(jù)作者們多年研究發(fā)現(xiàn),當(dāng)側(cè)體位于主體中部時(shí),不論是線性薄船理論還是Dawson型方法,都能得到與試驗(yàn)結(jié)果吻合較好的興波阻力曲線,而當(dāng)側(cè)體位于主體后部時(shí),線性薄船理論得到的結(jié)果與試驗(yàn)結(jié)果差別較大,因此選擇這種側(cè)體布局作為計(jì)算實(shí)例:該三體船側(cè)體縱向偏距為0.25 L,橫向偏距為0.2 L.

根據(jù)對(duì)稱性,半(右)側(cè)自由面計(jì)算區(qū)域縱向取為[-L,2L],橫向取為[0,L],即自由面網(wǎng)格船首方向伸展0.5 L,船尾方向伸展1.5 L,船寬方向伸展L.為對(duì)比分析,采用疊模流線網(wǎng)格和本文方法生成的自由面網(wǎng)格這兩類(lèi)自由面網(wǎng)格進(jìn)行計(jì)算,疊模流線網(wǎng)格即原始的Dawson型方法所采用的自由面網(wǎng)格,根據(jù)疊模計(jì)算得到的結(jié)果按流線劃分自由面網(wǎng)格,物理量的導(dǎo)數(shù)采用流線上的差分來(lái)表示.

流線網(wǎng)格取為均勻網(wǎng)格,L長(zhǎng)度內(nèi)劃分38個(gè)網(wǎng)格.由自由面疊模流線與平行于y軸的直線組成,如圖6所示,各流線在起點(diǎn)處y方向上間隔相等.三體船側(cè)體的存在會(huì)使得側(cè)體內(nèi)外流線 y方向上間隔不均勻,為此盡可能地保證同一寬度下各縱向位置處在y方向上的網(wǎng)格數(shù)接近相等,主體與側(cè)體之間橫向劃分網(wǎng)格數(shù)為 7.船體表面與自由面保持相同的網(wǎng)格縱向尺度,即主體縱向劃分 38個(gè)四邊形網(wǎng)格,側(cè)體縱向劃分 19個(gè)四邊形網(wǎng)格,主體和側(cè)體沿吃水方向上分別劃分10個(gè)和8個(gè)網(wǎng)格.

圖6 流線自由面網(wǎng)格Fig.6 Stream line free-surfacemesh

采用本文自由面網(wǎng)格快速生成系統(tǒng)生成兩種不同自由面上的“貼水線”網(wǎng)格進(jìn)行計(jì)算,兩者區(qū)別為傾角不同,一個(gè)傾角為 0°,即垂直網(wǎng)格(橫向網(wǎng)格線垂直于x軸線),類(lèi)似圖6的流線網(wǎng)格;另一個(gè)橫向網(wǎng)格線傾角為 45°,類(lèi)似圖 1所示網(wǎng)格,其他參數(shù)保持相同,即縱向單位船長(zhǎng)網(wǎng)格數(shù)為 38,橫向總網(wǎng)格數(shù)為 38,主體與側(cè)體間橫向網(wǎng)格數(shù)為 7,橫向網(wǎng)格密度控制因子為 1.0.如此參數(shù)設(shè)置是盡可能與流線網(wǎng)格保持一致,以客觀比較.

圖 7給出的是興波阻力系數(shù)計(jì)算結(jié)果與試驗(yàn)結(jié)果[8]的比較.圖中流線網(wǎng)格曲線是指采用流線網(wǎng)格興波阻力程序[9]根據(jù)圖 6所示流線網(wǎng)格計(jì)算得到的結(jié)果;0°網(wǎng)格曲線、45°網(wǎng)格曲線是指采用貼體網(wǎng)格興波阻力計(jì)算程序按上述貼體網(wǎng)格計(jì)算得到的結(jié)果,縱向?qū)?shù)采用迎風(fēng)格式計(jì)算,橫向?qū)?shù)采用中心差分格式計(jì)算,導(dǎo)數(shù)計(jì)算采用物理平面與計(jì)算平面的轉(zhuǎn)換關(guān)系在計(jì)算平面上進(jìn)行計(jì)算,假設(shè)計(jì)算平面坐標(biāo)(ξ,η)與物理平面坐標(biāo)(x,y)之間為一一對(duì)應(yīng)關(guān)系:

將物理空間的導(dǎo)數(shù)變換為對(duì)計(jì)算空間 ξ、η的導(dǎo)數(shù).

其中:

在計(jì)算平面上對(duì)自由面上物理量的導(dǎo)數(shù)進(jìn)行計(jì)算是本文三體船興波阻力數(shù)值計(jì)算方法的第 3個(gè)特點(diǎn),如前所述,前 2個(gè)特點(diǎn)是自由面網(wǎng)格生成技術(shù)方面的特點(diǎn).此外,圖 7還給出了根據(jù)線性薄船理論興波阻力計(jì)算結(jié)果[10-11]加以比較.

圖 7 興波阻力系數(shù)計(jì)算結(jié)果對(duì)比Fig.7 Wave drag comparison between the calculated results and experimental resu lts

從圖 7可以看出,與線性薄船理論計(jì)算結(jié)果相比,線性自由面條件的三維Rankine源方法與試驗(yàn)結(jié)果更加接近,薄船理論結(jié)果與模型試驗(yàn)結(jié)果偏差較大,由于算例側(cè)體位于主體后部,如前文所述,側(cè)體位于主體后部,伴有較強(qiáng)的噴濺應(yīng)該是計(jì)算結(jié)果與模型試驗(yàn)結(jié)果偏差的原因之一.

流線網(wǎng)格與 0°自由面網(wǎng)格計(jì)算結(jié)果差別甚微,這相互驗(yàn)證了流線網(wǎng)格、貼體網(wǎng)格興波阻力計(jì)算結(jié)果的正確性,同時(shí)也說(shuō)明采用自由面上的“貼水線”網(wǎng)格計(jì)算方法可以得到與流向網(wǎng)格方法一致的結(jié)果.流線網(wǎng)格生成費(fèi)時(shí)費(fèi)力,尤其是三體船興波阻力計(jì)算過(guò)程中自由面網(wǎng)格生成更加復(fù)雜,由于側(cè)體的存在,側(cè)體兩側(cè)緊貼側(cè)體表面的流線需要反復(fù)計(jì)算和調(diào)整才能得到,采用貼體網(wǎng)格預(yù)報(bào)三體船興波阻力計(jì)算過(guò)程則較為簡(jiǎn)潔方便.

45°自由面網(wǎng)格計(jì)算結(jié)果與流線網(wǎng)格、0°自由面網(wǎng)格計(jì)算結(jié)果有一定差別,與試驗(yàn)結(jié)果比較來(lái)看似乎 45°自由面網(wǎng)格計(jì)算結(jié)果與試驗(yàn)結(jié)果更加接近,分析可能原因是 45°自由面網(wǎng)格更符合興波的物理現(xiàn)象,但是否說(shuō)明斜網(wǎng)格一定優(yōu)于直網(wǎng)格尚待進(jìn)一步研究.

4 結(jié)束語(yǔ)

作者成功開(kāi)發(fā)了適用于三體船興波阻力計(jì)算的自由面網(wǎng)格快速生成系統(tǒng),該系統(tǒng)僅需簡(jiǎn)單輸入必要的船體型值和幾個(gè)控制參數(shù)即可快速生成三體船興波阻力計(jì)算的自由面貼體網(wǎng)格,并實(shí)現(xiàn)網(wǎng)格密度變化控制,同時(shí)可應(yīng)用于方尾三體船,避免了商用網(wǎng)格生成軟件生成三體船興波阻力計(jì)算自由面網(wǎng)格的繁瑣過(guò)程.對(duì)一艘Wigley數(shù)學(xué)船型興波阻力的計(jì)算說(shuō)明,采用系統(tǒng)生成的貼體網(wǎng)格計(jì)算三體船興波阻力,在計(jì)算過(guò)程簡(jiǎn)潔方便的同時(shí)還可以得出良好的計(jì)算結(jié)果.

自由面網(wǎng)格傾角變化對(duì)興波阻力計(jì)算結(jié)果有一定影響,需要進(jìn)一步研究.另外,方尾三體船興波阻力的計(jì)算與非方尾船略有不同,除需要該文已實(shí)現(xiàn)的生成方尾后流動(dòng)區(qū)域網(wǎng)格分塊外,還需在阻力計(jì)算程序中對(duì)方尾做相應(yīng)的特殊處理,這將在后續(xù)研究中考慮.

[1]高 斌,繆國(guó)平.關(guān)于Dawson方法中網(wǎng)格劃分的研究[J].海洋工程,2000,18(1):20-28.

GAO Bin,MIAO Guoping.Investigation of grid generation in Dawson method[J].Ocean Engineering,2000,18(1): 20-28.

[2]SAUNDERS JR P.An investigation of the resistance properties of a modern trimaran combatant ship base on Taylor standard series and series64[D].CA:naval Postgraduate School,1995:1-9.

[3]盧曉平,潘雨村.高速三體船興波阻力與片體布局優(yōu)化研究[J].水動(dòng)力學(xué)研究與進(jìn)展(A輯),2004,19(3): 347-359.

LU Xiaoping,PAN Yucun.An investigation of wave resistance on the high speed trimarans and their piece hull position optim ization[J].Journal of Hydrodynamics(Ser A),2004,19(3):347-359.

[4]WANG Hu,ZOU Zaojian.Numerical research on wavemaking resistance of trimaran[J].J Shanghai Jiao Tong University,2008,13(3):348-351.

[5]劉應(yīng)中.船舶興波阻力理論[M].北京:國(guó)防工業(yè)出版社,2003:89-90.

[6]陳京普,朱德祥,何術(shù)龍.雙體船/三體船興波阻力數(shù)值預(yù)報(bào)方法研究[J].船舶力學(xué),2006,10(2):23-29.

CHEN Jingpu,ZHU Dexiang,HEShulong.Research on numerical prediction method for wavemaking resistanceof catamaran/trimaran[J].Journal of Ship Mechanics,2006,10 (2):23-29.

[7]施法中.計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理 B樣條[M].北京:高等教育出版社,2001:1-11.

[8]BATTISTIN D,DANIELLIA,ZOTTI I.Numericaland experimental investigations on wave resistanceof trimaran congurations[C]//Proceedings of the 9th International Maritime Association ofMediterranean.Ischia,Italy,2000.

[9]DAWSON CW.A practical computer method for solving ship wave problems[C]//Proc of 2nd Intl Con f on Num Ship Hydrodyn.Berkeley,Calif,1977.

[10]王 中,盧曉平,鐘士崗.基于單元柯欽函數(shù)精確積分的多體船興波阻力計(jì)算[J].哈爾濱工程大學(xué)學(xué)報(bào), 2009,30(6):602-606.

WANGZhong,LU Xiaoping,ZHONG Shigang.Multi-hull ship wave making resistance calculation based on Kochin function integral analysis expression on the surface panels [J].Journal of Harbin Engineering University,2009,30 (6):602-606.

[11]PENG Hongxuan.Numerical computation of multi-hull ship resistance and motion[D].Halifax Dalhousie University,2001:9-38.