低傅氏數(shù)興波數(shù)值預(yù)報(bào)的不穩(wěn)定現(xiàn)象及其解決方案

梁 川，趙 峰

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

0 引 言

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

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

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

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

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

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

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

1 現(xiàn)象描述

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

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

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

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

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

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

圖3 常規(guī)設(shè)置Fr=0.138和Fr=0.28工況總阻力系數(shù)隨時(shí)間變化對(duì)比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）工況的計(jì)算出現(xiàn)了非常明顯的不穩(wěn)定現(xiàn)象。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 Fr=0.138工況不同設(shè)置下總阻力系數(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ù)隨時(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ù)增大。此時(shí)，單獨(dú)靠一種方式將很難控制計(jì)算穩(wěn)定性，需要多種方式聯(lián)合使用，例如在使用時(shí)間一階隱式的同時(shí)，還需要減小時(shí)間步長(zhǎng)。

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

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

4 結(jié) 論

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

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

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

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

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

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

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