帶自由液面有限長圓柱繞流數值模擬

陳松濤，趙偉文，萬德成*,，高洋洋

1 上海交通大學 船海計算水動力學研究中心，上海 200240

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

3 浙江大學 海洋學院，浙江 舟山 316021

0 引 言

隨著海洋油氣資源開發逐漸走向深遠海，各類深海平臺成為作業的主要裝備。圓柱形結構作為典型的結構之一，被廣泛應用于Spar平臺、TLP平臺以及半潛式平臺等深海結構物。在實際工程問題中，自由液面及自由端的存在極大地改變了結構物周圍的流動特性，從而對結構強度等造成額外的影響。相較于單相圓柱繞流問題[1-2]，自由液面的引入會在圓柱前駐點位置形成液面抬升，在尾流區域內形成液面凹陷，在遠場形成與船舶興波問題中類似的Kelvin波系。因此，自由液面的變形、翻卷及破碎會在一定程度上改變流場的流動特性。相較于無限長圓柱繞流問題，在自由端附近產生的漩渦將會影響典型的卡門渦街結構，帶來更為復雜的三維效應[3-4]。

近年來，國內外許多學者采用模型試驗和數值模擬的方法對該復雜問題展開了分析研究。Chaplin和Teigen[5]在拖曳水池中對不同吃水下的圓柱以雷諾數Re與弗勞德數Fr之間的固定比值（Re/Fr=2.79×105）開展了模型試驗，著重針對阻力系數和自由液面抬升進行了測量對比，發現圓柱前方的液面抬升與Re/Fr的取值范圍有關。Potts等[6]通過拖曳試驗對不同吃水比圓柱的水動力載荷和尾流場進行了分析，并采用3種不同的自由端形式進行了對比。Zhao等[7]在Fr=0.3，Re=4.2×104工況下對帶自由液面的無限長圓柱進行了數值模擬，并利用無因次化Rortex/Liutex渦識別方法從湍流結構的角度出發分析了自由液面給流動帶來的影響，結果表明在自由液面附近的尾流區域，原本交替脫落的渦街結構受到了較強的抑制作用。Koo等[8]采用大渦模擬（large eddy simulation，LES）方法和CLSVOF（coupled level set and volume of fluid）界面捕捉方法分析了不同Fr和Re數下圓柱繞流的自由液面效應，并通過與單相流動的對比，從流動分離和自由液面形態等角度展開了討論，發現Fr的影響更為明顯，且在高Fr下在圓柱前方觀察到了劇烈的波浪破碎現象。Rosetti等[9]采用非定常雷諾平均（unsteady Reynolds-averaged Navier-Stokes，URANS）方 法，在Fr=0.31，Re=4.3×104的工況下，對帶自由液面的低長徑比圓柱進行數值模擬，并通過與實驗的對比研究了自由液面及自由端的效應。Benitz等[10]針對Fr=0.65，Re=2 900工況下的有限長圓柱，利用不同的湍流模型進行數值模擬，主要分析了吃水比對受力和泄渦形式的影響。然而，先前研究所選取的工況均主要集中在低Fr和Re數。為了更好地研究其背后的規律和機理，需要借助數值模擬等手段對高Fr和Re數下的水動力特性進行分析與討論。

本文擬基于延時分離渦湍流模型和分段線性界面重構（PLIC）幾何類流體體積（VOF）法，利用自主開發的naoe-FOAM-SJTU求解器對Fr=1.1，Re=2.7×105工況下帶自由液面的有限長圓柱繞流開展數值模擬，并通過與實驗測量的自由液面結果進行對比來驗證求解器的可靠性。然后，在此基礎上對不同吃水位置處的受力、速度、流動分離和尾渦分布等進行對比分析，探討自由液面及自由端效應給流動帶來的影響，從而為實際的工程問題提供指導。

1 數值方法

本文采用的基于開源代碼框架OpenFOAM自主開發的求解器naoe-FOAM-SJTU[11]，其可靠性和精確度已在船舶阻力試驗[12]、敞水螺旋槳試驗[13]、船?槳?舵干擾[14]及海洋立管渦激振動（VIV）[15]/海洋平臺渦激運動（VIM）[16]等一系列船舶與海洋工程水動力學問題中得到廣泛驗證。

1.1 控制方程

對于兩相不可壓縮黏性流動，其控制方程包括連續性方程和動量方程，可寫成如下形式：

式中：U為速度；ρ為 加權密度；pd為動壓力；g為重力加速度；x為位置矢量；μeff為等效動力黏度；fσ為表面張力項；t為時間。

1.2 湍流模型

為了精確模擬高Re數下的大分離湍流流動問題，本文采用基于剪切應力輸運（shear-stress transport, SST）兩方程模型[17]的延遲分離渦模擬（delay detached-eddy simulation，DDES）技術。對于湍流流動，在工程實際應用中通常有2種方法：RANS方法和LES方法。RANS方法是將Navier-Stokes方程作時均化處理，將湍流中的物理量分解為時間平均項和脈動項，通過引入渦黏模型或雷諾應力模型來對方程組進行封閉，由于其濾去了時間脈動項，導致在求解大分離流動時往往很難獲得精細化的流場細節。LES方法是利用濾波函數將湍流運動分離成大尺度運動和小尺度運動，對于大尺度運動，直接進行求解，對于小尺度運動，則引入亞格子模型來模化其對大尺度運動的影響，通常需要在空間上保證足夠的網格分辨率以解析出高精度的湍流細節。本文采用的SSTDDES屬于一種混合雷諾平均/大渦模擬（hybrid RANS/LES）方法，即在近壁面的邊界層網格內采用RANS模式進行求解，在遠離壁面的分離流動區域則采用LES模式進行求解。相比LES，SST-DDES大大降低了近壁面附近的網格分辨率，同時也保證了尾流區域內的求解精度。其輸運方程包括湍流動能k和湍流耗散率ω方程，其形式如下：

式中：μ為介質的分子動力黏度；μt為渦流動力黏度；Pk為湍流動能產生項；F1為 一混合函數；CDkω為一限制函數；σk，σω， α，β 均為定常數；lDDES為DDES方法中的湍流長度，其定義如下：

式中：lRANS為 RANS模式下的湍流長度；CDES為一混合常數；fd為延遲函數；Δ為網格尺度。相較于原始的分離渦模擬（detached-eddy simulation，DES），DDES引入了延遲函數fd，可避免在近壁面處過早地從RANS模式切換到LES模式，定義如下：

式中：rd為 延遲因子；dw為到壁面的距離；S為應變率張量；Ω 為渦量；ν 為 介質的分子運動黏度；νt為渦流運動黏度；κ為von Karman常數；Cd1和Cd2為定常數。相關參數的具體定義和取值可參照文獻[18-19]。

1.3 界面捕捉方法

為了更好地捕捉該復雜流動中劇烈變化的自由液面，本文采用了OpenFOAM v8版本中最新植入的PLIC幾何類VOF方法[20]，其基本思想是：通過單個平面將網格單元切分成2部分，滿足每一時間步初始的相體積分數，并以此得到各個面的相體積分數通量，從而保證在求解相體積分數輸運方程時的精確性。圖1展示了該方法的二維示意圖，其主要流程為：

1） 判斷網格相體積分數αv， 若滿足 0<αv<1，則進行切割；

2） 將需要進行切割的網格相體積分數αv通過各頂點與網格中心處距離的倒數，插值得到各頂點的體積分數αp；

3） 根據各頂點體積分數αp的情況，判斷平面切割的位置，得到線性切割平面，如圖1中黑色虛線所示；

經過十年再就業工作的總結和思考，我們清醒地認識到，協解人員雖然和油田解除了勞動關系，但是從歷史淵源講、從感情的角度講、從利益紐帶的角度講，協解人員的石油情結沒有變。他們的生活、工作的主陣地仍在油田，他們始終都是“勝利人”，對他們的管理是油田躲不開、繞不過也推不出去的。做好協解人員的再就業和穩定工作是構建和諧油田的重要保障，是各級管理組織必須擔負的一項政治任務，任何時候都不能放松。那么，從心理學角度看，通過深入分析協解人員的思想和行為規律，積極運用心理干預技術進行有效嘗試，將有助于建立完善協解人員再就業工作長效穩定機制。

圖1 分段線性界面重構方法二維示意圖Fig.1 Two-dimensional diagram of PLIC method

4） 將切割得到的浸沒部分體積①除以網格體積②，得到相體積分數；

相比OpenFOAM以往采用的基于界面壓縮的代數類VOF方法，該方法提高了通量的計算精度，且同時能夠避免當壓縮系數設定過大時界面失真的問題。相比傳統的PLIC方法，其在拓展應用于三維非結構化復雜網格時，編程實現的復雜程度降低，且同時能夠與代數類方法進行結合，在一定程度上保證了求解復雜問題時的魯棒性，適用于一般的工程實際問題。

2 計算模型

本文參照美國愛荷華大學水利學研究所（IIHR）于2015年開展的拖曳水池實驗[21]，選取長徑比L/D=6（其中L為圓柱長，D為圓柱直徑）的圓柱作為研究對象。為了更好地體現自由液面對流動的影響，本文從中選取了Fr=1.1，Re=2.7×105的工況進行數值模擬。其中，均勻來流速度U∞參照實驗中的拖車速度設置為1.54 m/s。在此基礎上，水的密度為1 000 kg/m3, 運動黏度為1.14×10?6m2/s；空氣的密度為1 kg/m3, 運動黏度為1.48×10?5m2/s。

本文采用的計算域大小為40D×30D×10D，如圖2所示，其中圓柱置于水平面中心，吃水h=4D。圖3展示了本文采用的計算網格，由OpenFOAM自帶的snappyHexMesh工具繪制而成，總數量為8.3×106。其中，圖3（a）展示的是沿x-z平面方向的整體布局，圖3（b）中的紅色框線代表x-y平面上的局部加密區域，圖3（c）進一步展示了近圓柱表面的分級加密情況。圓柱表面共布置了8層邊界層網格，離壁面最近一層的網格高度滿足y+ <1。同時，為了更好地探究自由液面及自由端對流動產生的影響，在相應關注的區域進行了局部加密。

圖2 計算域Fig.2 Computational domain

圖3 計算網格Fig.3 Computational grids

上游入口處采用均勻來流速度入口邊界條件，下游出口處采用壓力出口邊界條件。圓柱表面設置為無滑移邊界條件，底部設置為對稱邊界條件，兩側及頂部均設置為零梯度邊界條件。時間步長為2×10?3D/U∞，滿足最大庫朗數Co始終小于1的要求。數值模擬的總時間為tU∞/D=120，取穩定的后半部分進行數據分析。

3 結果與分析

3.1 自由液面抬升

相比傳統的單相圓柱繞流問題，自由液面的引入使流動變為氣?液兩相流動。在一定來流速度下，自由液面會在圓柱前方產生抬升，在后方產生凹陷。圖4所示為不同縱剖面位置處的平均自由液面抬升結果，圖5所示為平均自由液面的等高線圖（圖中y= 0分割線以上的部分為實驗數據，以下部分為本文數值模擬結果）。由圖4和圖5可知，數值模擬結果與實驗值[21]吻合良好，對于圓柱前方的液面抬升，峰值誤差在3%以內，能夠很好地捕捉自由液面的變化情況，驗證了求解器的可靠性。

圖4 平均自由液面抬升剖線圖Fig.4 Profiles of the time-averaged free surface elevation

圖5 平均自由液面等高線圖Fig.5 Contour lines of the time-averaged free surface elevation

3.2 升阻力系數

式（8）和式（9）給出了兩相流下圓柱升、阻力系數的公式。為了更好地反映自由液面和自由端對于圓柱升、阻力在吃水方向上的影響，在數值模擬過程中將圓柱沿垂向（z方向）均分為了n份，每一份的長度為Δh=h/n（不應小于圓柱表面網格的最小尺度），每一份的受力為ΔF，在式（8）和式（9）中替換掉原有的F和h后，即可得到不同吃水深度位置處的升、阻力系數。

圖6和圖7所示分別為圓柱的平均阻力系數和升力系數均方根隨吃水變化的曲線，圖8和圖9所示為3個不同吃水深度的升、阻力時歷曲線。由圖6可以看出，平均阻力曲線在自由液面（x=0）及自由端位置（x= ?4）附近存在2個峰值，這可結合流場中的壓力與速度進行解釋。圖10所示為圓柱表面的時均壓力系數云圖（圖中，Cp為時均壓力系數）。從圖10中可以看出，圓柱后方產生的液面落差使得壓差相比吃水較深的位置更為明顯。圖11所示為流場x-z平面時均流向速度分布云圖（圖中，Uˉx/U∞為無因次化時均流向速度）。自由端的存在使得端部的局部速度增大，且存在明顯的上洗趨勢，故而在一定垂向范圍內改變了圓柱表面的低壓區域分布。從圖7可以看出，升力系數均方根僅在自由端附近存在1個明顯的峰值，這可結合下文的尾渦分析來看，這是因為該工況下較為劇烈的液面變形使得自由液面附近的泄渦形式變得復雜，不再類似于單相流動中交替脫落的卡門渦街。而對自由端而言，其端部額外產生的泄渦則增強了該位置附近的升力波動。從瞬時角度出發，圖8和圖9首先定性地驗證了上述時均數據的準確性，同時還可以額外發現，隨著吃水位置逐漸靠近自由端，升阻力系數的波動幅值均得到了極大的增加，表明該區域的流動存在著較強的三維瞬時特性。

圖6 平均阻力系數曲線Fig.6 Profile of the time-averaged drag coefficient

圖7 升力系數均方根曲線Fig.7 Root-mean-square profile of the lift coefficient

圖8 不同深度阻力系數時歷曲線Fig.8 Time histories of drag coefficient at different depths

圖9 不同深度升力系數時歷曲線Fig.9 Time histories of lift coefficient at different depths

圖10 圓柱表面時均壓力系數圖Fig.10 Time-averaged pressure coefficients of cylinder surface

圖11 x-z平面時均流向速度圖Fig.11 Time-averaged streamwise velocity on x-z plane

3.3 速度分布

圖12給出了圓柱近尾流場中不同橫剖面位置處的無因次化時均流向（x方向）和橫向（y方向）速度分布。從圖中可以看出，與深吃水位置相比，自由液面的存在延緩了順流向的速度“恢復”，且在遠離圓柱后方的一定區域內仍然有著較大的速度“損失”范圍。同時，橫向速度分布顯示，在近自由液面附近，流體呈現出向外側運動的趨勢，與吃水較深位置處的正好相反。這說明自由液面的引入在一定程度上改變了其附近的流動特性，今后需要與深吃水位置區分開來單獨予以討論。

圖12 時均速度剖面曲線Fig.12 Profiles of the time-averaged velocity

3.4 流動分離

圖13所示為圓柱表面不同吃水位置處的時均速度梯度曲線，其中迎流方向的前駐點位置定義為0°。從圖中可以看出，圓柱表面的分離角θ隨吃水呈先減小后增大的趨勢。為了更好地分析該現象，圖14給出了圓柱表面不同吃水位置處的時均壓力系數曲線。從圖中可以看出，隨著吃水的增加，圓柱表面前半部分的正壓梯度逐漸增大。相應地，這對于靠近自由液面的吃水位置處（z=?1D），減小了其后半部分的逆壓梯度，從而延遲了流動分離的發生。對于靠近自由端的吃水位置處（z=?3D），雖然其前半部分的正壓梯度較大，但其背部的低壓區域較其余位置也更為顯著，使得壓力恢復減緩，這同樣導致了逆壓梯度的減小。

圖13 時均速度梯度曲線Fig.13 Profiles of the time-averaged velocity gradient

圖14 時均壓力系數曲線Fig.14 Time-averaged pressure coefficient profiles

3.5 尾渦分析

圖15所示為自由液面及不同吃水位置處的瞬時垂直渦量（ωz）云圖，選取的時刻為tU∞/D=120。從圖中可以看出，自由液面上的自由剪切層在圓柱后方呈向外擴張的形式，這與3.3節中提到的橫向速度分布規律相符合。同時，在自由液面上占主導的是由液面變化產生的細碎漩渦。隨著吃水的增加，自由剪切層不再主要受自由液面的影響，逐漸向內側收縮，且此時能夠在尾流場發現交替脫落的卡門渦街。但隨著吃水位置逐漸靠近自由端，垂直渦量不再沿順流向傳播，而是主要集中于圓柱背部。圖16所示為采用無因次化Rortex/Liutex渦識別方法[22]繪制的圓柱尾流場三維渦結構圖。相比其他傳統的渦識別方法，Rortex/Liutex渦識別方法能夠更好地排除剪切運動的“污染”，且能同時捕捉尾流場中的強弱渦。從圖16中可以看出，自由端產生的卷擰狀漩渦從端部脫落后斜向上運動，這在一定程度上抑制了卡門渦街沿順流向的發展。

圖15 瞬時垂直渦量云圖Fig.15 Instantaneous vertical vorticity

圖16 瞬時渦結構Fig.16 Instantaneous vortical structure

4 結 論

本文利用自主開發的naoe-FOAM-SJTU求解器，對帶自由液面的有限長圓柱繞流進行了數值模擬。在湍流模型上，選擇SST-DDES方法，在自由液面捕捉上，采用PLIC幾何類VOF方法。通過將自由液面抬升與實驗結果進行對比，驗證了求解器的可靠性，并對不同位置處的受力、速度、流動分離和尾渦分布進行了分析，主要得到如下結論：

1）自由液面和自由端的存在使得平均阻力系數沿吃水方向存在2個對應的峰值，其主要原因是圓柱前后液面落差形成的壓差和經端部加速后的上洗流動；升力系數均方根僅在自由端附近存在1個明顯的峰值，其主要原因是較為劇烈的液面變形抑制了泄渦的產生，而端部則產生了明顯的卷擰狀漩渦。

2） 自由液面的引入使得流向速度相比深吃水位置“恢復”延緩，并使橫向速度呈現出明顯向外側流動的分布趨勢，后面的瞬時渦量云圖則從側面驗證了這一點。

3） 圓柱表面的流動分離在自由液面和自由端的雙重作用下，均由于逆壓梯度的減小呈向后延遲的變化規律，前者與前半部分減小的正壓梯度有關，后者則與背部的低壓區域有關。

4） 基于無因次化Rortex/Liutex渦識別方法，自由液面上占主導的為由液面變化產生的細碎漩渦，而自由端產生的卷擰狀漩渦則在一定程度上抑制了卡門渦街沿流向的發展。

這說明在實際工程問題中，自由液面和自由端的存在會使流動特性發生較大變化，這就需要根據實際情況分析所帶來的影響。目前，本文僅對單個工況進行了數值模擬來驗證該計算模型的準確性，今后還需繼續補充不同工況下的計算，以進行系統化的對比研究。