基于GEKO 湍流模型的JBC 流場數值模擬

許輝，日野孝則，陳作鋼*

1 上海交通大學 海洋工程國家重點實驗室，上海 200240

2 上海交通大學 船舶海洋與建筑工程學院，上海 200240

3 橫濱國立大學 工學研究院，日本 橫濱 240-8501

0 引 言

在船舶水動力性能預報中，除經典的船模試驗和經驗方法外，計算流體力學（computational fluid dynamics，CFD）方法也得到了廣泛應用。在黏性作用和湍流效應不可忽略的流場中， N-S 方程是流體運動的控制方程。雷諾平均納維-斯托克斯（Reynolds-averaged Navier-Stokes， RANS）數值模擬方法是采用經驗性較強的湍流模型來對平均后的N-S 方程進行封閉，是目前工程計算中最常用、也最經濟的N-S 方程數值求解方法。

SSTk-ω 湍流模型具有良好的壁面距離自適應性和逆壓梯度流預報能力，但因其在構造過程中引入了經驗性參數，使得湍流模型不具有普適性。Menter[1]提出了GEKO（Generalizedk-ω）湍流模型的構造思路，該模型在SSTk-ω 湍流模型的基礎上，通過引入和調整不同的流動控制參數拓展了其適用范圍。

JBC（Japan bulk carrier）是一種設計裝備有導管型節能裝置的好望角型散貨船，被作為2015 東京CFD 研討會的測試船型。Hino 等[2]對JBC 船型進行拖曳水池試驗及風洞試驗，測定了船體阻力及伴流場數據，并給出了JBC 船型CFD 計算的相關建議。隨后針對2015 東京CFD 研討會提交的JBC 算例，Hino 等[3]又進行了系統、全面的分析，結果顯示各類CFD 方法在阻力預報方面表現均相對較好，不同的湍流封閉方法對伴流場的預報影響均較明顯。Sun 等[4]、Lidtke 等[5]和Schuiling等[6]均采用傳統的兩方程模型對包含和不包含節能導管的JBC 流場進行了數值模擬，結果顯示船體阻力預報數據與試驗結果吻合良好，伴流預報給出的舭部縱向渦強度偏小。 Korkmaz[7]采用各向異性的EASM 模型對包含和不包含節能導管的JBC 流場進行了數值模擬，結果顯示船體阻力預報結果略小于試驗數據，而伴流預報中的舭部縱向渦強度則略大于試驗數據，整體來說吻合較好。除了采用傳統的經典湍流模型外，也嘗試采用其他數值方法對JBC 流場進行了相關研究。Bensow 等[8]嘗試采用部分平均方法（partially averaged Navier-Stokes，PANS）對JBC 流場進行了數值預報，其通過采取調整控制參數的方式來控制RANS 方法的應用范圍，結果顯示隨著RANS 計算范圍的減小，船體的黏性阻力逐漸減小，整體的阻力預報也逐漸偏小，船尾伴流場的預報結果與單純的RANS 方法相近，認為混合方法中的近壁湍流場并沒有及時轉捩和發展。Queutey 等[9]采用分離渦模擬-剪切應力輸運（DES-SST）方法對JBC 模型的流場進行了瞬態計算，結果表明JBC 的舭部渦是由一些強非定常性的湍流相干結構疊加而成，DES 方法通過還原該類渦結構得到了與試驗結果一致性很好的船尾湍動能分布。隨后，Kornev 等[10]同樣采用混合雷諾平均-大渦模擬（RANS-LES）方法對JBC 流場進行了數值模擬，計算結果顯示采用混合方法預報的摩擦阻力Cf明顯偏小，而壓差阻力Cp明顯偏大，同時還研究分析了船體附近湍流結構的發展過程，得到了相比RANS 方法與試驗結果吻合更好的伴流場速度及湍動能分布預報數據。總體而言，JBC 船型作為一種標準算例已進行過大量研究，是驗證GEKO 湍流模型性能的良好算例。

本文將闡述GEKO 湍流模型的部分構造理論并采用JBC 船型進行計算驗證，以檢驗構造思路的可行性和新模型對工程CFD 計算的適用性，為湍流模型工程應用的改進提供參考。

1 RANS 數值模擬理論

1.1 控制方程

在船舶水動力性能數值預報中，在流體黏性和湍流作用不可忽略的條件下，流動現象控制方程主要采用三維不可壓縮黏性流動的RANS 方程組：

式中：ui，uj為瞬時速度分量；xi，xj（i，j=1，2，3）為3 個方向的空間坐標；ρ 為水的密度；p為壓力；ν 為運動黏性系數；fi為體積力；t為時間；有上劃線的物理量為平均值，有單引號的物理量為脈動值。方程右端第3 項為不封閉項，表現為2 個脈動速度分量的統計相關：

也即雷諾應力，是湍流模型理論的核心。

1.2 SST k-ω 湍流模型

渦黏湍流模型采用Boussinesq 渦黏假定[11]來表達雷諾應力，將脈動量的統計相關表達為平均量的函數：

式中：μT為動力渦黏系數，相應的運動渦黏系數νT=μT/ρ；δij為Kronecker 算子；k為單位流體質量湍動能；Sˉi j為平均應變張量。渦黏湍流模型作為線性各向同性湍流模型，構造簡單，計算中魯棒性強，在工程上得到了廣泛應用。

引入渦黏假定后，湍流模型理論的核心轉變為封閉渦黏系數μT。

Menter[12]在大量工程試驗和計算的基礎上提出了一種經驗方法，即將標準k-ε 模型[13]與Wilcoxk-ω 模型[14]進行混合：

式中：F1為標準k-ε 模型與Wilcoxk-ω 模型的混合函數； φ1為Wilcoxk-ω 模型中的任意封閉系數；φ2為標準k-ε 模型中的任意封閉系數； φ為混合后的任意封閉系數。通過封閉系數混合，可以得到BSL （Baseline）模型：

式中，β*，β，σk，γ，σω，σω2均為經驗性封閉系數；μ為流體動力黏性系數。

SSTk-ω 湍流模型在BSL 模型的基礎上重定義渦黏系數μT：

式中：a1為經驗系數；Ω為渦量大小；F2為修正剪切應力輸運效應的函數。式(6)～式(8)即為SSTk-ω模型的k方程、ω 方程和渦黏系數封閉式。相比Wilcoxk-ω 模型和標準k-ε 模型，SSTk-ω 模型在逆壓梯度流中有著良好的預報能力，具有壁面距離自適應性，在工程計算中運用非常廣泛。

1.3 GEKO 湍流模型

Menter[1]根據試驗和工程計算經驗提出，湍流模型封閉系數的細微改動會引起計算結果精確度的大幅度提升或是惡化，可以利用該思路為不同工程應用場景設計最適用的封閉系數。為便于進行系數調節，Menter[1]在提出GEKO 湍流模型時引入了多個調節參數，其中流動分離參數CS和旋轉修正參數CRC對船舶流場的影響較大。但Menter沒有公開其調節參數的構造方法，本文中CS及CRC采用的是Hino 等[15]研究中的定義。

首先，流動分離參數CS通過修正式(5)中的封閉系數 φ1來達到控制流動分離的目的。修正后的封閉系數 φ為

對比式(9)和式(5)可知，當CS= 1.75 時，φ=φ1，式(5)未進行任何修正，此時，GEKO 湍流模型等價于SSTk-ω 湍流模型；當CS= 1.00 時， φ=φ2=φ，此時，GEKO 模型的封閉系數全部為標準k-ε模型設置。在本文計算中，CS的取值在1.00～1.75之間。根據Menter[1]提供的參考算例，在其經過計算及試驗數據標定的范圍0.7＜CS＜ 2.5 內，隨著CS的增大，邊界層流動對逆壓梯度應該逐漸敏感，從而流動分離也隨之加劇。

GEKO 模型中，旋轉修正參數CRC的設計思路源于Hellsten[16]對原始SSTk-ω 模型的修正，目的是提高SSTk-ω 模型對流動旋轉和壁面曲率的敏感性。在GEKO 模型中，旋轉修正參數CRC是通過將式(7)中右側第2 項（即耗散項βρω2）乘以修正函數F4來實現對模型的控制，F4定義為

其中：

式中：Ri為用于曲面流動修正的Richardson 數；S為平均應變率；為平均渦量張量。式(10)中，旋轉修正參數CRC主要用來表達修正幅度，不同的CRC會對結果造成不同程度的影響[17]，Hellsten建議CRC取為3.6。本文的CRC取值在0～3.6 之間。通過引入CS和CRC，應用于船舶流動的GEKO 湍流模型即構造完成。不過，模型參數對不同流動工況的適用性還需要進行大量的研究來予以確定，本文只進行了初步的嘗試。

2 數值模擬

2.1 幾何模型

本文以JBC 船型為研究對象對肥大船型進行GEKO 湍流模型計算測試，模型示意圖如圖1所示。JBC 船型的主尺度如表1 所示。計算中，以光體模型為計算對象，數值計算和對比阻力試驗均在縮尺比為1∶40 的模型尺度進行，對應的試驗數據來自日本海上技術安全研究所（national maritime research institute，NMRI）的系列試驗[2]。

圖1 JBC 船體模型示意圖Fig. 1 Schematic diagram of the JBC bare hull model

表1 JBC 船型主尺度Table 1 Main dimensions of JBC ship

2.2 網格生成

計算網格采用O-O 型拓撲全結構網格，由商業軟件Pointwise生成，計算域范圍為-1.5 ≤x/Lpp≤ 4.5，-2.0 ≤y/Lpp≤ 0，-2.0 ≤z/Lpp≤ 0，如 圖2 所示。計算域網格單元總數約為142.5 萬。為解析計算近壁面處邊界層的分離等流動細節，第1 層網格節點被設置于y+～1 的空間位置，垂向網格增長率設置為1.1。船艉部分的網格分布如圖3 所示。

圖2 計算域及邊界條件示意圖Fig. 2 Schematic diagram of computational domain and boundary conditions

圖3 JBC 船艉部分網格示意圖Fig. 3 Grid near the stern of JBC ship

2.3 求解器及計算設置

數值計算采用SURF ver. 7 求解器[18]，該求解器為基于有限體積法的三維不可壓縮流動RANS 方程數值求解器，其采用二階迎風格式離散動量方程，同時采用人工可壓縮方法處理壓力速度耦合。計算中，來流的雷諾數與試驗中相同，取以垂線間長Lpp為特征長度的雷諾數，即Re=7.46 × 106。試驗中采用的弗勞德數Fr= 0.142。為集中考察GEKO 湍流模型的計算效果，自由面采用對稱條件，略去興波阻力的影響。

在計算過程中，流動分離參數CS和旋轉修正參數CRC各有3 種取值，所有的計算設置如表2所示。

表2 各工況對應參數設置Table 2 Parameter setting of varying cases

根據Menter 的計算標定，流動分離參數CS的取值可以在0.7 ＜CS＜ 2.5 范圍內選取。當CS= 1.75時，GEKO 湍流模型等價于SSTk-ω 湍流模型，即C3，C3R1，C3R2 是在SST 模型的基礎上進行的不同程度的CRC修正；當CS= 1.00 時，GEKO 模型封閉系數全部為標準k-ε 模型設置，即C1，C1R1，C1R2 是在標準k-ε 模型的基礎上進行的不同程度的CRC修正；當CS= 1.5 時，GEKO 模型為過渡模型，即其既不是SSTk-ω 模型，也不是標準k-ε模型。流動分離參數CS的取值可以超出1.00～1.75 范圍，以進一步抑制或者加強湍流模型對流動中逆壓梯度的敏感性，但從物理意義上來說，因其無既有模型可以對比，故本文暫不涉及。

對于CRC的取值，當CRC= 3.6 時，GEKO 湍流模型采用Hellsten[16]建議的旋轉效應修正幅度；當CRC= 0 時，GEKO 模 型 采 用 無 旋 轉 效 應 的 修正；當CRC= 2.75 時則為過渡值，用于考察不同修正幅度對計算結果的影響。

3 計算結果對比及分析討論

3.1 阻力系數對比

在阻力系數對比中，總阻力系數Ct由摩擦阻力系數Cf和黏壓阻力系數Cp這2 部分組成，并由數值模擬計算直接得出這2 部分的結果。試驗結果中的Cf由ITTC 1957 經驗公式給出，Ct由試驗[2]測量的形狀因子1+k=1.314 與Cf相乘獲得。試驗及所有的數值模擬結果如圖4 所示。

圖4 數值預報阻力分量與試驗數據對比Fig. 4 Comparisons between resistance components of numerical predictions and experimental data

1） 流動分離參數CS作用分析。

以C1，C2，C3 工況對比為例，當旋轉修正參數CRC= 0 時，隨著流動分離參數CS的增大，船體總阻力系數Ct也逐漸增大，其中黏壓阻力系數Cp略有增大，摩擦阻力系數Cf則是顯著增大。當CRC為其他取值時，Ct及其分量隨CS的變化趨勢相同。

從理論上分析，當船舶表面的邊界層流動分離加劇時，黏壓阻力應有增大的趨勢，這與計算結果所表現的基本吻合，證明流動分離參數CS在黏壓阻力調節方面基本合理。對于摩擦阻力部分，JBC 船型的摩擦阻力部分是隨CS的增大而增大的。

在其他相關研究[15]中，KVLCC2 船型的摩擦阻力系數則是隨CS的增大而減小，JBC 與KVLCC2同為肥大型船型，卻表現出了與CS不同的關系。研究中所用湍流模型的理論構建和程序開發均是基于Menter 對其所設計GEKO 湍流模型的性能描述，因目前的測試算例相對有限，研究中CS所用的線性插值可能難以滿足Menter 對所有模型的要求。最終，由于JBC 與KVLCC2 型線的區別，使得在構造較簡單的湍流模型時表現出了不同的性能。

2） 旋轉修正參數CRC作用分析。

以C3，C3R1，C3R2工況對比為例，此時CS= 1.75，GEKO 模型等價于進行不同CRC修正的SSTk-ω模型。當流動分離參數CS保持恒定時，隨著旋轉修正參數CRC的增大，總阻力系數Ct幾乎沒有變化，其中摩擦阻力部分略有減小，黏壓阻力部分略有增大。當CS為其他取值時，總阻力系數Ct及其分量也有相同的變化趨勢。總體而言，當CS恒定時，各種CRC設置下的阻力預報結果變化不大。

在其他研究[15]中，對于阻力系數與CRC的關系， KVLCC2 船型也表現出了與JBC 船型相同的變化趨勢。從理論的角度分析，對于肥大型船，湍流模型針對流動旋轉效應和船體幾何曲率的敏感性所產生的阻力影響可能很有限。

3.2 船尾伴流場及壓力分布對比

GEKO 模型的參數設置會對流場細節造成一定的影響。對于JBC 船型，槳盤面位于x/Lpp= 0.986 4處，在試驗中對該處的伴流場進行了測量[2]，如圖5所示。各參數設置下的數值模擬結果如圖6 所示。肥大型船槳盤面的典型特征是由縱向渦引起的伴流扭曲，這在所有計算結果中均有體現。湍流模型的選取是決定船尾伴流場的關鍵因素，采用渦黏假定的各向同性線性湍流模型預報的縱向渦鉤狀低速區通常小于試驗測定值[3]。不同設置下GEKO 模型預報的船尾附近壓力分布及伴流場分布如圖7 所示。圖中，船體表面為無量綱壓力系數P的分布。

圖5 試驗中槳盤面處（x/Lpp=0.986 4）伴流場速度分布Fig. 5 Velocity distribution of wake field at the propeller plane(x/Lpp=0.986 4)

圖6 數值預報的槳盤面(x/Lpp = 0.986 4)處伴流場速度分布對比Fig. 6 Comparisons of velocity distribution of wake fields at the propeller plane (x/Lpp = 0.986 4) by numerical prediction

圖7 各參數設置下船艉壓力系數及伴流場對比Fig. 7 Comparisons of the pressure distributions and wake fields near the stern at varying parameter settings

式中：p∞為來流壓力；U∞為來流速度；g為重力加速度。船尾附近為中縱剖面的無量綱速度U的分布，其中U=u/U∞，u為流場的縱向速度分量。

1） 流動分離參數CS作用分析。

以C1，C2，C3 工況對比為例，當CRC= 0 時，隨著CS的增大，圖6 所示槳盤面伴流場的鉤狀低速區明顯增大，但在各種設置下均未捕獲到試驗中U= 0.2 的低速區。總體而言，當CS= 1.75 時，數值計算結果與試驗數據的吻合程度最高。相對應的，圖7 中縱剖面上的低速區隨著CS的增大也有增大趨勢；由船體壓力分布來看，船體后肩部舭部低壓區的面積是逐漸增大的，艉部高壓區的變化則相對較小，舭部附近的等壓線隨著CS的增大逐漸匯集，因而從舭部到船艉的逆壓梯度隨著CS的增大具有增大的趨勢，從而加劇了流動分離。當CRC為其他取值時，數值預報的結果也有類似的變化趨勢。

從理論上分析，隨著CS的增大，GEKO 模型對逆壓梯度更加敏感，流動分離則應更加劇烈，因而必然會對伴流場產生明顯的影響。而船體后肩處低壓區的擴散很可能是因為流動分離的加劇削弱了船艉壓力的恢復，從而導致黏壓阻力增大。

在其他研究[15]中，KVLCC2 船型隨CS也有類似的變化趨勢，該船型與試驗結果吻合最好的槳盤面伴流場出現在CS=1.85 時，超出了本文CS的取值范圍。經對比可知，盡管當CS＞1.75 時沒有既有的湍流模型可以對比，但仍有進一步研究的必要。

2） 旋轉修正參數CRC作用分析。

以C3，C3R1，C3R2 工況對比為例，此時，GEKO模型在SSTk-ω 模型的基礎上進行的CRC修正。當CS= 1.75 時，隨著CRC的增大，圖5 中槳盤面伴流場的鉤狀低速區逐漸增大并產生輕微的上移和向中縱剖面的靠攏。總體而言，當CRC= 3.6 時，GEKO 模型預報的槳盤面伴流場與試驗數據的吻合程度最高。而圖6 所示縱剖面上的伴流場低速區隨CRC的變化則不明顯，主要集中在槳帽附近；由船體壓力分布來看，船體舭部低壓區隨著CRC的增大有面積增大的趨勢，同時艉部的高壓區幾乎沒有變化，舭部附近的等壓線隨著CRC的增大逐漸匯集，而其他區域的變化則不明顯。由于船體舭部的曲率較大，在該區域的流場細節隨CRC變化明顯的同時，其他區域幾乎沒有變化，表明CRC的增大改變了模型對幾何曲率的敏感性。當CS為其他取值時，其數值預報結果也有類似的變化趨勢。對比來看，盡管CRC對JBC 船型阻力預報結果的影響有限，但對伴流場預報結果的影響則相對明顯。

在其他研究[15]中，KVLCC2 船型槳盤面伴流場的低速區域范圍隨CRC的變化不明顯，主要表現為低速區的輕微上移。從理論上分析，隨著湍流模型對流動旋轉效應和船體幾何曲率敏感性的增大，很可能會導致肥大型船舶的縱向渦強度增大，進而導致船體舭部區域的低壓區域擴散和黏壓阻力的輕微上升。

4 結 論

湍流模型理論中包含有大量的量綱分析過程和經驗性因素，在很大程度上還沒有觸及到湍流的物理機理本質，因而不具有普適性。然而，在數值計算中若根據實際工況調整湍流模型計算參數的設定，采用經濟的RANS 方法仍有可能使計算獲得相對合理、精確的預報結果，而采用GEKO湍流模型通過設定新參數來對傳統湍流模型的封閉系數進行系統性調節是一種快捷的方法。本文嘗試將GEKO 模型應用到了JBC 船型的流場數值預報。文中所構造的GEKO 湍流模型引入了流動分離參數CS和旋轉修正參數CRC，用于調節流動分離強度和流動旋轉效應。CS和CRC對數值計算的影響需要通過大量的數值試驗進行測試，本文僅針對JBC 船型的流場預報進行了初步的嘗試，主要結論如下：

1） 對于JBC 船型的阻力預報，CS的增大會導致摩擦阻力的顯著增大和黏壓阻力的小幅增大；CRC對船體阻力的影響較小，CRC的增大會導致摩擦阻力的略微減小和黏壓阻力的略微增大，總阻力值基本保持不變。

2） 對于JBC 船型的伴流場預報，CS和CRC的影響均比較明顯。隨著CS的增大，槳盤面伴流場的鉤狀低速區會產生相應程度的增大，中縱剖面的低速區也會產生相應程度的增大；隨著CRC的增大，槳盤面低速區在增大的同時會產生輕微的上移和向中縱剖面的靠攏，而中縱剖面的伴流場變化不明顯。

3） 對于JBC 船型的船艉壓力預報，CS的增大會導致船體后肩低壓區面積的增大，進而削弱船艉的壓力恢復并增大黏壓阻力，同時，還會增大由舭部至船艉的逆壓梯度，加劇流動分離；CRC的增大則主要會導致幾何曲率較大舭部區域的流場細節變化，其他區域的變化則不明顯。

4） GEKO 模型的不同參數設置對數值計算結果具有直接影響，可在一定程度上提高數值預報適用性和精確度的潛力。

在GEKO 湍流模型新參數的合理選取和對不同計算工況的適應性方面，今后仍需采用大量的數值計算進行探索。