低傅氏數興波數值預報的不穩定現象及其解決方案

梁 川，趙 峰

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

0 引 言

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

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

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

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

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

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

本文將針對低傅氏數工況興波數值預報不穩定的現象展開相關研究，嘗試分析其背后的數值機理，并提出相應的解決方案。

1 現象描述

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

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

表1 阻力計算常規參數設置Tab.1 Regular CFD settings for resistance computation

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

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

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

圖3為Fr=0.138和Fr=0.28工況總阻力系數隨時間變化曲線。可以看到，Fr=0.28時波動的幅值更小且呈逐漸縮小趨勢，而Fr=0.138時波動振幅很大（通常Froude數越低，振幅越小），且波峰波谷的包絡線呈震蕩變化。

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

上述種種現象表明，在同樣采用常規設置情況下，低傅氏數（Fr=0.138）工況的計算出現了非常明顯的不穩定現象。

2 控制方程離散格式數值分析

出現上述現象是因為在比較低的傅氏數下計算，其計算穩定性不如高傅氏數情況，而計算穩定性與控制方程的離散格式的數值特性密切相關，特別是差分格式的數值耗散特性，具體可以參見相關教材[10-12]。

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

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

針對控制方程（式（2）），按照常規設置的離散格式（時間二階隱式格式、對流二階迎風格式、擴散項二階中心差分格式、源項采用顯式處理）進行數值耗散項（離散格式和原微分形式之間差異中的偶數階微分項）分析。

以時間項為例，二階隱式差分格式與原微分項的差值為

暫不關注其他項，不難推導出式（2）的數值耗散形式，如式（6）所示：

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

3 源于理論指導的數值實踐

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

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

減小逆耗散：從逆耗散項的表達式中可以看到，減小時間步長Δt可以減小逆耗散。

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

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

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

另外，從正常的低傅氏數流動的波形范圍和波高幅值來看，傅氏數越低，波形范圍越小（波長越短，且主要集中在船首附近）、波高幅值也越小。因此，在進行網格劃分時，為了能夠捕捉到自由面處的流動，自由面位置常需要布置更加細密的網格，具體為水平方向重點對船體首尾區域加密，同時垂直方向也需進行適當加密。由于本文常規設置的自由面網格已經采用了較好的加密布置，本文未進行自由面網格加密的數值對比實踐。

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

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

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

另外，由自由面產生的逆耗散項形式可知，當Froude 數持續降低，逆耗散項將持續平方指數增大。此時，單獨靠一種方式將很難控制計算穩定性，需要多種方式聯合使用，例如在使用時間一階隱式的同時，還需要減小時間步長。

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

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

4 結 論

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

（1）對于阻力計算這類穩態問題，適合于采用時間一階隱式格式；

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

（3）適當加密自由面附近的網格；

（4）綜合應用上面三種策略。

需要說明的是，由于控制策略的具體細化與所采用的求解器有很強的關聯性，當遇到低傅氏數計算不穩定問題時，本文所給的策略將起到一個方向指引的作用。

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