梁 川，趙 峰

（中國船舶科學(xué)研究中心，江蘇無錫214082）

0 引 言

船舶航行于水面，并伴隨有相應(yīng)的船行波，自由液面的存在是船舶性能數(shù)值預(yù)報的一個顯著特征。很長一段時間內(nèi)，船舶粘性流動和自由面的計算是分開考慮的，即用勢流理論方法處理自由面興波，而另外通過求解RANS 方程來計算船體的粘性邊界層，這種分離處理的方法，不能考慮粘性和自由面之間的耦合作用。隨著計算機(jī)性能的快速提升和相關(guān)CFD 算法的深入發(fā)展，特別是諸如VOF、Level Set 等自由面捕捉方法的發(fā)展和成熟，使得同時考慮粘性和自由液面影響的船舶CFD 技術(shù)獲得了長足的進(jìn)步，取得了眾多應(yīng)用成果[1-7]，并愈來愈朝著數(shù)值水池應(yīng)用愿景方向發(fā)展[8-9]。

傅氏（Froude）數(shù)（式（1））是一個衡量流體慣性力與約化重力（流體質(zhì)點(diǎn)受到重力與浮力的合力，在其他位置為零，在自由面存在一個階躍）比值的無量綱參數(shù)，是一個區(qū)分船舶性能與類別的重要參數(shù)，其區(qū)分定性作用類似于航空中的馬赫(Mach)數(shù)。

在采用CFD 技術(shù)進(jìn)行船舶水動力性能研究時，其難點(diǎn)主要在Froude 數(shù)范圍的兩端，即高傅氏數(shù)（Fr＞0.5）和低傅氏數(shù)（Fr＜0.15）工況。

高傅氏數(shù)船舶常見于一些高速快艇，其計算難點(diǎn)主要在于船體興波作用很大、船體姿態(tài)變化劇烈、船體運(yùn)動與興波作用強(qiáng)烈耦合等非線性特征突出，相關(guān)的CFD技術(shù)在不斷研究和發(fā)展中。低傅氏數(shù)船舶常見于大型油輪和超大散貨船，隨著船舶綠色環(huán)保理念的提升，運(yùn)輸船舶越來越向大型化方向發(fā)展，對于低傅氏數(shù)船舶性能的數(shù)值預(yù)報也越來越具有顯著的工程價值。低傅氏數(shù)工況的計算難點(diǎn)主要在于隨著Froude 數(shù)越來越低，計算將變得越來越不穩(wěn)定，且這一不穩(wěn)定現(xiàn)象，在自研求解器或者商用求解器（如CFX，F(xiàn)luent和StarCCM+等）中普遍存在。

在低傅氏數(shù)船舶粘流數(shù)值模擬方面，已有不少學(xué)者進(jìn)行了相關(guān)研究。典型的如Toxopeus 等[5]，統(tǒng)計研究了四款主流求解器（MARIN 水池的ReFRESCO、Iowa大學(xué)的CFDSHIP、ECN 大學(xué)的ISIS-CFD 和商用軟件STAR-CCM+）對標(biāo)模KVLCC2 做偏航運(yùn)動模擬（聚焦于Fr=0.064 2 工況）的流動情況，其中，ReFRESCO 和CFDSHIP 未考慮自由面影響，ISIS-CFD 和STAR-CCM+對自由面進(jìn)行了考慮，不過文中未提及考慮自由面后可能出現(xiàn)的計算不穩(wěn)定現(xiàn)象。

目前，關(guān)于該不穩(wěn)定現(xiàn)象的文獻(xiàn)報導(dǎo)很少，分析和解釋這一現(xiàn)象的更少。這主要是當(dāng)Froude數(shù)較低時，自由面興波的作用較小，通常會當(dāng)做疊模（自由面處設(shè)置為對稱邊界條件）來進(jìn)行處理。然而，針對該問題的深化研究除了理論上解釋關(guān)于某個計算現(xiàn)象的疑問外，同時也有其實(shí)際應(yīng)用價值：（1）低傅氏數(shù)工況，興波阻力雖較小但并不表示不存在，而且自由面的存在對粘壓阻力存在一定影響，比如水池試驗(yàn)獲取船舶形狀因子的方法主要在低傅氏數(shù)段進(jìn)行，目前ITTC 推薦方法經(jīng)常失效的一些特殊情況（球艏貼近水面，方尾浸沒），均與自由面對粘壓阻力的影響有關(guān)；（2）對于船舶在波浪中低速航行的工況，自由面將無法采用疊模處理。

本文將針對低傅氏數(shù)工況興波數(shù)值預(yù)報不穩(wěn)定的現(xiàn)象展開相關(guān)研究，嘗試分析其背后的數(shù)值機(jī)理，并提出相應(yīng)的解決方案。

1 現(xiàn)象描述

筆者在進(jìn)行DTMB5415標(biāo)模（CFD WorkShop 2010標(biāo)模[4]，高、中、低三個典型航速算例分別對應(yīng)Fr=0.41、0.28 和0.138）阻力計算時，發(fā)現(xiàn)低Froude 數(shù)工況計算反常，留心記錄并與相關(guān)業(yè)內(nèi)人士交流討論后，發(fā)現(xiàn)這個現(xiàn)象在船舶CFD 模擬中普遍存在，是一類共性問題。其主要表現(xiàn)為：（1）計算波形十分反常，且隨計算時間增加興波范圍和幅度越來越大；（2）計算的阻力曲線波動很大，總體來看是計算不穩(wěn)定的典型特征。

由于本文所關(guān)注的低傅氏數(shù)工況興波數(shù)值預(yù)報不穩(wěn)定現(xiàn)象在諸多CFD求解器中普遍存在，因此，關(guān)于本文所采用的CFD 方法只作一個簡要介紹：基于求解不可壓RANS方程的CFD 求解器、計算網(wǎng)格生成采用重疊網(wǎng)格技術(shù)、自由面處理采用Level Set方法。在數(shù)值實(shí)踐初期，采用平時初步總結(jié)出的一套參數(shù)設(shè)置進(jìn)行計算，也稱為常規(guī)設(shè)置，如表1所示。下面以標(biāo)模DTMB5415在Fr=0.138和Fr=0.28時舉例說明相關(guān)現(xiàn)象。

表1 阻力計算常規(guī)參數(shù)設(shè)置Tab.1 Regular CFD settings for resistance computation

圖1 Fr=0.138工況常規(guī)設(shè)置下不同時間步的波形對比Fig.1 Comparison of wave contours at different time steps with the regular settings when Fr=0.138

圖2 Fr=0.28工況常規(guī)設(shè)置下不同時間步的波形對比Fig.2 Comparison of wave contours at different time steps with the regular settings when Fr=0.28

圖1 為Fr=0.138 工況，常規(guī)設(shè)置下不同時間步的波形對比圖，隨著時間的推移，波動的范圍和幅度越來越大，波面的能量越來越積聚，這與真實(shí)的低傅氏數(shù)興波很小的特征明顯不符；圖2 為Fr=0.28 工況的波形對比圖，可以看出，波形變化穩(wěn)定，與物理直觀相符。

圖3為Fr=0.138和Fr=0.28工況總阻力系數(shù)隨時間變化曲線?？梢钥吹剑現(xiàn)r=0.28時波動的幅值更小且呈逐漸縮小趨勢，而Fr=0.138時波動振幅很大（通常Froude數(shù)越低，振幅越小），且波峰波谷的包絡(luò)線呈震蕩變化。

圖3 常規(guī)設(shè)置Fr=0.138和Fr=0.28工況總阻力系數(shù)隨時間變化對比Fig.3 Comparison of Ct curve over time steps with the regular settings when Fr=0.138 and Fr=0.28

上述種種現(xiàn)象表明，在同樣采用常規(guī)設(shè)置情況下，低傅氏數(shù)（Fr=0.138）工況的計算出現(xiàn)了非常明顯的不穩(wěn)定現(xiàn)象。

2 控制方程離散格式數(shù)值分析

出現(xiàn)上述現(xiàn)象是因?yàn)樵诒容^低的傅氏數(shù)下計算，其計算穩(wěn)定性不如高傅氏數(shù)情況，而計算穩(wěn)定性與控制方程的離散格式的數(shù)值特性密切相關(guān)，特別是差分格式的數(shù)值耗散特性，具體可以參見相關(guān)教材[10-12]。

需要通過離散差分來近似表達(dá)各微分項(xiàng)，離散格式和原微分形式之間的差異（更高階的數(shù)值耗散或者數(shù)值色散項(xiàng)）將影響數(shù)值計算的特性。特別是偶數(shù)階的數(shù)值耗散項(xiàng)，將很大程度上影響計算的穩(wěn)定性。

一般來說，采用隱式格式或者迎風(fēng)格式將產(chǎn)生數(shù)值耗散（數(shù)值粘性），增加計算穩(wěn)定性；采用顯式格式則會產(chǎn)生數(shù)值逆耗散，增加計算不穩(wěn)定性，如控制方程中不同變量的分離求解（即通過采用上一時間步的結(jié)果，假設(shè)其他變量已知，然后輪流求解當(dāng)前變量，并通過不斷迭代達(dá)到收斂），將引入不穩(wěn)定因素。在不可壓N-S 方程求解中，一般采用分離求解的策略，如式（2）中的源項(xiàng)，即采用顯式處理，從而引入不穩(wěn)定因素。

針對控制方程（式（2）），按照常規(guī)設(shè)置的離散格式（時間二階隱式格式、對流二階迎風(fēng)格式、擴(kuò)散項(xiàng)二階中心差分格式、源項(xiàng)采用顯式處理）進(jìn)行數(shù)值耗散項(xiàng)（離散格式和原微分形式之間差異中的偶數(shù)階微分項(xiàng)）分析。

以時間項(xiàng)為例，二階隱式差分格式與原微分項(xiàng)的差值為

暫不關(guān)注其他項(xiàng)，不難推導(dǎo)出式（2）的數(shù)值耗散形式，如式（6）所示：

另外，同理易得時間項(xiàng)采用一階隱式的數(shù)值耗散形式，如式（7）所示：

3 源于理論指導(dǎo)的數(shù)值實(shí)踐

與Froude 數(shù)相關(guān)的逆耗散項(xiàng)有使計算發(fā)散的趨勢，因此數(shù)值實(shí)踐的思路可圍繞如何來抑制這一趨勢展開，包括兩個方面，一是增加數(shù)值正耗散，另一個是減小數(shù)值逆耗散。

增加正耗散：一個是增加時間項(xiàng)耗散，比如時間差分由二階隱式改為一階隱式，實(shí)際上對于阻力計算這種穩(wěn)態(tài)問題，用一階格式比較合適；另一個是增加對流項(xiàng)耗散，比如由二階格式改為一階格式，但是這種方法對計算結(jié)果影響較大，要慎用。

減小逆耗散：從逆耗散項(xiàng)的表達(dá)式中可以看到，減小時間步長Δt可以減小逆耗散。

圖4 為在常規(guī)設(shè)置基礎(chǔ)上，分別將時間項(xiàng)改為一階隱式（圖4（c））、將時間步長縮小5 倍取Δt=0.002（圖4（b））與常規(guī)設(shè)置（圖4（a））在相同來流時間下的波形對比圖，另將三者在y=0.08Lpp處的波形曲線繪制于圖4（d）?？梢钥吹剑统Ｒ?guī)設(shè)置相比，大范圍的不合理波動均得到了不同程度的抑制，其中，采用時間一階隱式的方案在船體周圍展現(xiàn)的合理波形更加豐富，縮小5倍時間步長的方案在船體外圍仍存在一些不合理的波動。仔細(xì)分析數(shù)值耗散公式（式（6））會發(fā)現(xiàn)：減小時間步長既會減小逆耗散，同時也會減小時間項(xiàng)正耗散；單純增加正耗散（時間二階隱式改一階隱式）效果最好。

圖4 Fr=0.138工況不同參數(shù)設(shè)置的波形對比圖Fig.4 Comparison of wave contours with the different settings when Fr=0.138

圖5 是Fr=0.138 工況不同設(shè)置下總阻力系數(shù)隨時間變化曲線，可以看到：與常規(guī)方案相比，采用時間一階隱式的方案，波動的幅度大幅縮小，并最終趨于穩(wěn)定；從總阻力系數(shù)的平均值來看，兩者相差無幾，且與試驗(yàn)值吻合良好。

另外，從正常的低傅氏數(shù)流動的波形范圍和波高幅值來看，傅氏數(shù)越低，波形范圍越?。úㄩL越短，且主要集中在船首附近）、波高幅值也越小。因此，在進(jìn)行網(wǎng)格劃分時，為了能夠捕捉到自由面處的流動，自由面位置常需要布置更加細(xì)密的網(wǎng)格，具體為水平方向重點(diǎn)對船體首尾區(qū)域加密，同時垂直方向也需進(jìn)行適當(dāng)加密。由于本文常規(guī)設(shè)置的自由面網(wǎng)格已經(jīng)采用了較好的加密布置，本文未進(jìn)行自由面網(wǎng)格加密的數(shù)值對比實(shí)踐。

圖6 是Fr=0.28 工況，時間一階隱式和常規(guī)設(shè)置下總阻力系數(shù)隨時間變化曲線，可以看到曲線的波動幅度相差不大，一階隱式略小，這與Fr=0.138 工況差異明顯，從側(cè)面驗(yàn)證了計算穩(wěn)定性與Froude數(shù)平方呈反比的理論推導(dǎo)，也指出低傅氏數(shù)的不穩(wěn)定現(xiàn)象尤其值得關(guān)注。

圖5 Fr=0.138工況不同設(shè)置下總阻力系數(shù)隨時間變化曲線Fig.5 Comparison of Ct curves over time steps with the different settings when Fr=0.138

圖6 Fr=0.28工況不同設(shè)置下總阻力系數(shù)隨時間變化曲線Fig.6 Comparison of Ct curves over time steps with the different settings when Fr=0.28

另外，由自由面產(chǎn)生的逆耗散項(xiàng)形式可知，當(dāng)Froude 數(shù)持續(xù)降低，逆耗散項(xiàng)將持續(xù)平方指數(shù)增大。此時，單獨(dú)靠一種方式將很難控制計算穩(wěn)定性，需要多種方式聯(lián)合使用，例如在使用時間一階隱式的同時，還需要減小時間步長。

圖7 為在Fr=0.1 時，時間均采用一階隱式、時間步長分別為Δt=0.01 和Δt=0.005 時，總阻力系數(shù)隨時間變化曲線。可以看到：在Δt=0.01時，即使采用時間一階隱式格式，總阻力曲線也會出現(xiàn)比較大幅度的高頻振蕩；當(dāng)Δt=0.005 時，曲線才變得光滑。當(dāng)Froude 數(shù)繼續(xù)降低，經(jīng)測試，F(xiàn)r=0.05 時，在采用一階隱式時間格式的同時，時間步長需取Δt=0.002才能獲得光滑的結(jié)果。

圖7 Fr=0.1工況不同設(shè)置下總阻力系數(shù)隨時間變化曲線Fig.7 Comparison of Ct curves over time steps with the different settings when Fr=0.1

4 結(jié) 論

本文針對船舶在低傅氏數(shù)興波數(shù)值預(yù)報中容易出現(xiàn)計算不穩(wěn)定的普遍現(xiàn)象展開了相關(guān)研究，從控制方程離散差分格式出發(fā)，推導(dǎo)出船舶CFD所特有的（帶自由面）不穩(wěn)定項(xiàng)形式，從理論上來解釋這一現(xiàn)象。同時，根據(jù)所推導(dǎo)的數(shù)值耗散項(xiàng)的表達(dá)形式，進(jìn)行了充分的數(shù)值實(shí)踐研究，提出了低Froude數(shù)計算穩(wěn)定性的控制策略供相關(guān)從業(yè)人員參考借鑒：

（1）對于阻力計算這類穩(wěn)態(tài)問題，適合于采用時間一階隱式格式；

（2）減小時間步長能夠有效抑制發(fā)散；

（3）適當(dāng)加密自由面附近的網(wǎng)格；

（4）綜合應(yīng)用上面三種策略。

需要說明的是，由于控制策略的具體細(xì)化與所采用的求解器有很強(qiáng)的關(guān)聯(lián)性，當(dāng)遇到低傅氏數(shù)計算不穩(wěn)定問題時，本文所給的策略將起到一個方向指引的作用。

致謝：本文的完成也得益于與吳乘勝研究員、李勝忠研究員、邱耿耀高級工程師、蔣一工程師以及其他相關(guān)科研人員的啟發(fā)性交流，作者在此對他們的貢獻(xiàn)表示由衷的感謝！