黃衛(wèi)剛，邱遼原，姜治芳，陸 超，李 鵬

中國艦船研究設計中心，湖北武漢 430064

0 引 言

船舶繞流場數(shù)值預報是當前船舶計算流體力學研究的一個重點，該技術正發(fā)展成一種具備較高精度、可用于相似船舶航行性能比較的技術手段。采用CFD技術預報船舶流場存在的難點是：船舶運動時的流體動力會導致船舶出現(xiàn)升沉以及縱傾等運動；高質(zhì)量的網(wǎng)格劃分以及網(wǎng)格差異對計算的影響；數(shù)值預報與模型試驗之間的誤差問題。

船舶運動時出現(xiàn)的升沉以及縱傾等變化是導致數(shù)值預報結果與模型試驗結果之間出現(xiàn)誤差的主要原因之一，且這種誤差隨著航速的增加而越發(fā)明顯，如何處理兩者之間的差異成為數(shù)值預報中需首先解決的問題。國內(nèi)外對此已開展了各種研究，提出了諸如采用動網(wǎng)格等技術來解決該問題的措施。但在實際工程應用中，動網(wǎng)格技術對計算條件的要求較高。一般情況下，要采用動網(wǎng)格開展船舶繞流數(shù)值模擬，需并行幾十個CPU開展計算，整個計算周期較長，計算過程中的網(wǎng)格變形有可能會導致網(wǎng)格質(zhì)量變差，尤其是當計算模型的網(wǎng)格越復雜時，這種情況出現(xiàn)的可能性就越大，而網(wǎng)格質(zhì)量變差有可能會影響預報精度。如何在只有少量CPU且不采用動網(wǎng)格等技術的前提下實現(xiàn)較高精度的船舶繞流數(shù)值預報是本文要解決的主要問題。

數(shù)值模擬時的環(huán)境條件：計算機為4核CPU、內(nèi)存為8 G、數(shù)值模擬軟件的授權數(shù)為4個。本文的研究結論均是以此計算限制條件為基礎而開展。首先，以KCS和KVLCC為對象標定計算模型，然后以DTMB 5415為對象開展網(wǎng)格類型及規(guī)模對計算的影響分析，提出在計算條件受限的條件下宜采用六面體網(wǎng)格開展船舶繞流數(shù)值預報的建議。為進一步提高預報精度，本文提出了一種預先給定升沉的船舶繞流數(shù)值預報技術，并將以某條船為對象對該技術進行驗證性計算。最后，以兩條船為對象開展數(shù)值分析，進一步驗證CFD技術在船型優(yōu)化中的效用。

1 國內(nèi)外研究情況

模型尺度的船舶繞流模擬正在不斷進步，隨著技術的進步，船舶計算對象考慮的也將越來越真實。

在自由液面模擬方面，自由液面跟蹤方法非常有效。Gorski［1］以一條穿浪內(nèi)傾船為對象，采用自由液面捕捉的方式進行了數(shù)值模擬，數(shù)值模擬確定的沿船體表面及船尾的興波與模型試驗測定的興波之間吻合較好。Raven等［2］說明了自由液面捕捉技術仍有提升的空間，通過消除方程中的時間項并嚴格解決穩(wěn)態(tài)自由液面粘性流的方式，可以通過較少的數(shù)值計算來確定興波，該興波與系列60以及Dyne水池的模型試驗結果之間的吻合程度較高。在2007年的國際數(shù)值船舶水動力學會議上，西班牙Polytechnic大學的研究者指出，采用六面體網(wǎng)格計算船舶繞流能取得較高的精度。

當傅汝德數(shù)較低時，因船舶運動引起的升沉以及縱傾較小，因而可以直接采用不計及升沉及縱傾的方法來計算船舶阻力，此時，模型試驗的阻力與數(shù)值預報的阻力之間誤差較小。當傅汝德數(shù)增加時，船舶運動產(chǎn)生的升沉以及縱傾逐漸明顯，此時，若還采用不計及升沉及縱傾的方法計算阻力，則該阻力與模型試驗的阻力間的誤差就會逐漸明顯，在Fn較高時可能會引起約10%甚至更高的差異。因此，在當前階段，許多船舶粘性流數(shù)值模擬仍采用預先設定升沉和縱傾的方式。例如，Azcueta［3］闡述了采用VOF方法可以高精度預報一條單體瀕海戰(zhàn)斗艦在 Fn=0.25～0.9時的阻力情況；Azcueta［4］證明了一條高速艇在 Fn=0.75～4.0間的相似結論。

2 計算原理

1）應用RANS方程和連續(xù)性方程解決湍流問題時，由于增加了雷諾應力項而使基本方程不封閉，因而需要根據(jù)湍流的運動規(guī)律構建合適的湍流模型以解決方程的封閉性問題。為了更好地處理高應變率及流線彎曲度比較大的流動，本項目湍流模型采用了SST模型。

2）由于自由液面需要考慮空氣和水的相互作用，因而屬多相流問題。本研究處理自由液面問題時使用的是VOF方法。這種方法是一種以流體占據(jù)網(wǎng)格單元體積份額的途徑來跟蹤自由面演化的方法，可以跟蹤發(fā)生復雜變形的自由面。因VOF方法的流體質(zhì)量守恒在方程的構建中就得到了滿足，容易擴展到三維，體積分數(shù)的演變只需要用到相鄰的單元，因此得到廣泛應用。

3 計算條件受限情況下排水型中、高速船流場預報存在的問題及解決措施

3.1 存在的問題

在計算條件受限的情況下，排水型中、高速船流場預報存在的主要問題有：

1）船舶運動時的流體動力導致船舶出現(xiàn)升沉以及縱傾等運動。

排水型高速船與低速船在航行過程中存在明顯差異。排水型高速船在運動過程中的流體動力較明顯，導致船舶在運動過程中存在較大的縱傾和升沉，而低速船的升沉和縱傾則相對較小。很顯然，在固定吃水以及升沉的前提下，船舶在運動過程中形成的升沉和縱傾會成為影響數(shù)值預報結果與模型試驗結果間差異的主要因素，且會隨著傅汝德數(shù)的增加而更加明顯。

表1所示為一條縮尺比為1∶22的高速集裝箱船船模隨著傅汝德數(shù)的增加其升沉和縱傾的變化情況。在實船狀態(tài)，當Fn=0.21時，該船的升沉約為80 mm，在 Fn=0.39時，升沉約為360 mm?？梢婋S著傅汝德數(shù)的增加，升沉和縱傾的變化也會比較明顯，并將成為影響數(shù)值預報結果與模型試驗結果間誤差的關鍵因素。

表1 某船模在不同傅汝德數(shù)下的升沉和縱傾變化Tab.1 The sinkage and trim of a ship at different Fn

2）高質(zhì)量的網(wǎng)格劃分以及網(wǎng)格差異對計算的影響。

3）數(shù)值預報與模型試驗之間的誤差問題。

3.2 解決措施研究

為了解決在計算條件受限情況下排水型中、高速船流場預報存在的問題，本文提出了下列解決措施：

1）在軟、硬件條件受限的情況下，建議直接采用六面體網(wǎng)格，且網(wǎng)格質(zhì)量應滿足Angle＞20°，Aspect ratio＜100，Maximum Mesh Expansion Factor≤20等要求。該措施能保證在計算條件受限時數(shù)值預報結果具有較高的精度，減少對計算條件的要求。

2）如果設計船舶的傅汝德數(shù)較小，可以直接采用不計及升沉和縱傾的方法來計算裸船體阻力。如果設計船舶的傅汝德數(shù)較大，則建議采用預先給定升沉的方法，該措施能在避免使用動網(wǎng)格等技術的前提下解決固定吃水時船模所受浮力和重力不平衡的問題，有效提高計算精度。

為了驗證上述措施的可行性，本文以DTMB5415，KCS和KVLCC等船模為對象開展了數(shù)值預報，DTMB5415，KCS及KVLCC這三型標模的尺度以及三維模型如表2和圖1～圖3所示。

表2 標模尺度簡表Tab.2 Standard model’s dimension

圖1 DTMB 5415模型Fig.1 The model of DTMB5415

圖2 KCS模型Fig.2 The model of KCS

圖3 KVLCC模型Fig.3 The model of KVLCC

3.2.1 計算模型標定

計算模型滿足下列條件：

1）網(wǎng)格為六面體網(wǎng)格，Angle＞20°，Aspect ratio ＜100，Maximum Mesh Expansion Factor≤20。

2）湍流模型為SST模型。

3）流體域范圍：為避免遠方邊界條件對近船體流場的干擾，在文獻［5-10］的基礎上，計算區(qū)域入口取為船體艏部向上游延伸至2倍船長處；出口取為艉部向下游延伸至2倍船長處；區(qū)域外邊界分別為由對稱面（船體縱中剖面）向右舷方向和設計水線面向下延伸約0.8倍船長處；區(qū)域上邊界取為設計水線面向上約2倍吃水高度處。

4）邊界條件：入口為Inlet；出口為Opening；主船體壁面為No Slip Wall；流體域的底面、頂面及側面為Free Slip Wall；流體域的中線面為SYMMETRY。

在滿足上述條件的前提下，對標模KCS和KVLCC開展了數(shù)值預報。表3和表4所示為數(shù)值預報的裸船體阻力與模型試驗的裸船體阻力間的誤差對比。

表3 KCS模型試驗和數(shù)值預報的阻力對比Tab.3 The resistance’s difference between KCS’s model test and CFD

表4 KVLCC模型試驗和數(shù)值預報的阻力對比Tab.4 The resistance’s difference between KVLCC’s model test and CFD

表3和表4表明，在限制條件下，采用六面體網(wǎng)格預報的裸船體阻力與模型試驗阻力之間的誤差可以控制在工程精度范圍內(nèi)。這驗證了各種文獻所報道的結論，即在航速較低時，采用六面體網(wǎng)格預報的流場與模型試驗之間的吻合程度較好。

3.2.2 網(wǎng)格類型及規(guī)模對計算的影響

為進一步分析網(wǎng)格差異對計算的影響，針對DTMB 5415模型，分別對流體域為簡單四面體網(wǎng)格、四面體網(wǎng)格+邊界層網(wǎng)格以及六面體網(wǎng)格的單流體域開展了數(shù)值預報。受內(nèi)存所限，六面體網(wǎng)格的網(wǎng)格數(shù)量為240萬，四面體網(wǎng)格的網(wǎng)格數(shù)量約為450萬，四面體網(wǎng)格+邊界層網(wǎng)格的網(wǎng)格數(shù)量約為480萬。計算的網(wǎng)格如圖4所示。

圖4 DTMB 5415網(wǎng)格劃分情況Fig.4 The mesh of DTMB 5415

計算速度Vm=2.096 m/s，該狀態(tài)下船模的模型試驗阻力為45.1 N。不同的網(wǎng)格對CFD計算裸船體阻力誤差的影響的計算結果如表5所示。

表5 不同網(wǎng)格對數(shù)值預報阻力的誤差影響Tab.5 Influence of different meshes on the ship resistance calculation results

由表5可知，六面體網(wǎng)格的阻力預報結果與模型試驗結果之間的誤差最小，四面體網(wǎng)格+邊界層網(wǎng)格的次之，四面體網(wǎng)格的最大。

為進一步分析造成這種差異的原因，對尾部螺旋槳位置的x方向伴流進行了對比，如圖5所示。

由圖5可知，在計算條件受限的前提下，與模型試驗測得的沿主船體的流場形態(tài)相比，六面體網(wǎng)格計算確定的流場形態(tài)與模型試驗結果最接近，四面體+邊界層網(wǎng)格計算確定的流場形態(tài)次之，四面體網(wǎng)格計算確定的流場形態(tài)與模型試驗測得的流場形態(tài)間的誤差最大。出現(xiàn)該現(xiàn)象的原因為：

圖5 網(wǎng)格差異對尾部螺旋槳位置x方向伴流影響Fig.5 Influence of different meshes on the x-direction flow around ship’s propeller

1）在進行數(shù)值繞流的邊界層計算時，六面體網(wǎng)格沿船體表面分布的網(wǎng)格厚度可以很小，Y+可以小于200甚至是接近于10，能滿足湍流模型的需求，有助于精確捕捉邊界層的流場。四面體+邊界層網(wǎng)格相對于六面體網(wǎng)格而言，其沿船體表面分布的網(wǎng)格能滿足湍流模型的需求，但在邊界層與流體域之間的連接區(qū)域，其網(wǎng)格變化較大，導致離邊界層較遠區(qū)域的流場形態(tài)與模型試驗間的差異較大。四面體網(wǎng)格相對于六面體網(wǎng)格而言，其沿船體表面分布的網(wǎng)格厚度最大，離滿足湍流模型的需求尚存在一定的差距。

2）在開展自由液面計算時，六面體網(wǎng)格在自由液面附近能設置較多的網(wǎng)格，可為捕捉自由液面提供良好的條件。四面體+邊界層網(wǎng)格雖然彌補了四面體網(wǎng)格在捕捉邊界層附近流場上的缺點，但由于自由液面附近的網(wǎng)格較少，影響了自由液面計算的精度，從而導致計算的自由液面與模型試驗確定的自由液面之間的差異較大。四面體網(wǎng)格的自由液面捕捉情況與四面體+邊界層網(wǎng)格的情況類似。

3.2.3 預先給定升沉的船舶繞流數(shù)值預報技術

如前文所述，傅汝德數(shù)的增加會導致數(shù)值預報結果與模型試驗結果間的誤差逐漸加大，因此，本文接下來將要解決減少這種誤差的方法。一般情況下，可以采用動網(wǎng)格等技術來解決，但這種方法存在計算時間較長、計算過程中容易出現(xiàn)發(fā)散等問題，在計算條件受限的情況下實施難度很大。針對這種情況，這里將介紹一種預先給定升沉的船舶繞流數(shù)值預報方法。

根據(jù)計算排水量狀態(tài)下對應的吃水來計算船舶阻力。在正常情況下，計算1000～2000步時，計算的垂向力（即船舶所受浮力）就會出現(xiàn)如圖6所示的情況。這時，可以根據(jù)船體所受垂向力與模型重量間的差值來重新確定吃水，然后，根據(jù)新的吃水重新開展計算，直至船體垂向力與模型重量間的差異小于0.5%為止。此時，就可以認為數(shù)值預報時的排水量狀態(tài)與模型試驗的排水量狀態(tài)是一致的，一般迭代2～3次就可以得到計入升沉影響的阻力。

現(xiàn)以某型集裝箱船為例分析計入升沉與不計入升沉對數(shù)值預報結果的影響。該船模型重量為14400 N，當Fn=0.28時，如果采用不計入升沉的船舶繞流數(shù)值方法開展計算，此時計算確定的浮力為14120 N，比模型重量小了約280 N。如果采用計入升沉的船舶繞流數(shù)值方法開展計算，計入升沉后的浮力為14400 N，與模型的重量一致。當計入升沉的影響后，船舶的吃水增加，所受的浮力也隨之增加，從而達到重量與浮力的平衡，因此，此時的升沉狀態(tài)便與模型試驗狀態(tài)基本一致。

同時，還對比了計入升沉與不計入升沉對阻力的影響。當Fn=0.28時，根據(jù)模型試驗確定的阻力為120 N。計入升沉與不計入升沉后的阻力計算結果如表6所示。

表6 計入升沉與不計入升沉對阻力計算的影響Tab.6 The resistance calculation results with or without sinkage taken into account

由表6可知，不計入升沉時的阻力相對模型試驗確定的阻力120 N小了約3.7%，計入升沉后的阻力相對模型試驗的阻力小了約1.0%。計入升沉后的排水量更接近于模型試驗狀態(tài)，此時兩者之間的阻力誤差也更小。計入升沉的影響能有效減少數(shù)值預報結果與模型試驗結果之間的誤差，尤其是當船舶的傅汝德數(shù)較大時，這種效果更加明顯。當兩型存在較大差異的船進行阻力對比時，建議采用文中所述的方法。

4 驗證性計算

研究CFD技術的最終目的是應用到船型優(yōu)化中。為了進一步檢驗現(xiàn)有計算經(jīng)驗的可靠性，在限制條件下，針對兩個集裝箱船的相似方案（方案A和方案B）進行了流場計算。方案A和方案B均為長度約8 m的相似模型，船型系數(shù)和主尺度基本一致。計算的網(wǎng)格為六面體網(wǎng)格，采用的湍流模型為SST模型。數(shù)值預報的阻力對比如表7和表8所示。

由表7、表8可知：

1）在方案A和方案B中，流場計算確定的阻力與模型試驗之間的阻力誤差均小于5%。

2）兩型船的數(shù)值預報趨勢與模型試驗趨勢基本一致。

3）數(shù)值預報阻力與模型試驗阻力的差異符合下列趨勢：隨著傅汝德數(shù)的增加，模型試驗阻力與數(shù)值預報阻力的誤差逐漸加大。

表7 方案A計算阻力和試驗結果比較Tab.7 The project A resistance’s difference between model test and CFD

表8 方案B計算阻力和試驗結果比較Tab.8 The project B resistance’s difference between model test and CFD

5 結 論

本文探討了計算條件限制下船舶繞流數(shù)值預報精度的解決措施，認為要提高數(shù)值預報精度可以采取以下措施：

1）在開展船舶繞流數(shù)值預報時，如果其軟、硬件環(huán)境與本文類似，直接采用六面體網(wǎng)格比較合適。

2）如果所設計船舶的傅汝德數(shù)較小，可以直接采用固定吃水和縱傾的方法來計算裸船體阻力。一般情況下，數(shù)值預報結果與模型試驗結果之間的誤差能控制在工程允許的精度范圍內(nèi)。

3）如果所設計船舶的傅汝德數(shù)較大，此時船舶的升沉和縱傾已成為影響數(shù)值預報與模型試驗間差異的主要因素時，建議采用文中提出的預先給定升沉的方法。這種方法雖然會增加一些計算時間，但在計算條件等限制情況下，相對動網(wǎng)格等技術而言更簡單、實用。

通過本文的驗證性計算還表明，在相同的網(wǎng)格參數(shù)以及邊界條件設置下，采用數(shù)值預報來判別型線變化對阻力以及流場形態(tài)影響的方法是可行的。

［1］GORSKI J J.Present state of numerical ship hydrodynamics and validation experiments［J］.Journal of Offshore Mechanics and Arctic Engineering，2002，124：74-80.

［2］RAVEN H C，PLOEG A P，STARKE A R.Computation of free-surface viscous flows at model and full scale by a steady iterative approach［C］//Proceeding 25th Symposium Naval Hydrodynamics.Canada，2004.

［3］AZCUETA R.Steady and unsteady RANSE simulations for planning crafts［C］//The 7th International Conference on Fast Sea Transportation，F(xiàn)AST2003.Italy，2003.

［4］AZCUETA R.Steady and unsteady RANSE simulations for littoral combat ships［C］//Proceeding 25th Symposium on Naval Hydrodynamics.Newfoundland and Labrador，Canada，2004.

［5］王中，盧曉平.基于細分曲面的興波阻力計算方法［J］.海軍工程大學學報，2006，18（增刊）：20-25.WANG Zhong，LU Xiaoping.A calculation method for ship wave making resistance based on subdivision surface theory［J］.Journal of Naval University of Engineering，2006，18（Supplement）：20-25.

［6］PENG Q J，HOU Y J.A VOF based numerical model for breaking waves in surf zone［J］.Chinese Journal of Oceanology and Limnology，2006，24（1）：57-64.

［7］GAO Q X.Numerical simulation of free surface flow around ship hull［J］.Journal of Ship Mechanics，2002，6（3）：1-13.

［8］ZHANG Zhirong，ZHAO Feng，LI Baiqi.Numerical calculation of viscous free-surface flow about ship hull［J］.Journal of Ship Mechanics，2002，6（6）：10-16.

［9］HINO T.A study of grid dependence in navier-stokes solutions for free surface flows around a ship hull［J］.Journal of the Society of Navy Architects of Japan，1994，176：11-18.

［10］SHIOTANI S，KODAMA Y.Numerical computation of viscous flows with free surface around a series 60 model［J］.Journal of the Society of Naval Architects of Japan，1996，180：39-50.