基于IDDES方法的沖刷地形海底管道繞流數值模擬

曲孟祥，楊 博，梁丙臣, 2，劉 欣，張 嶔

(1.中國海洋大學 工程學院，山東 青島 266100；2.山東省海洋工程重點實驗室，山東 青島 266100)

海底管道對于近海油氣工業中碳氫化合物的運輸十分重要。在可侵蝕泥沙床面上，波浪和水流的作用會導致管道附近的泥沙沖刷，由于受到流動、地形以及土質等各種環境因素的不同影響，沖刷現象十分復雜。隨著沖刷坑的發展，當流體繞過管道時，在管道周圍會形成復雜且不穩定的分離和尾流結構，并且在管道后方產生旋渦脫落現象，引起管道的渦激振動，最終導致管道疲勞損壞，從而造成環境污染和重大經濟損失[1]。因此，研究沖刷過程中圓管周圍水動力特征對于石油和近海工業具有重要意義。截止到目前，對于海底管道繞流的研究主要分為兩個方面，一是對沖坑地形圓管周圍水動力的研究，二是對靠近平坦海床的圓管周圍水動力的分析。

對沖坑地形管道周圍流場的研究最早可追溯到20世紀70年代，Mao[2]通過試驗研究了管道和沖刷床面之間的相互作用，觀察了在穩定水流及波動條件下不同雷諾數、Shields參數情況下圓管周圍沖刷情況，并得到了沖刷初始階段5個特征階段的地形剖面。Jensen等[3]根據Mao[2]的試驗結果，對沖刷過程中的5個特征階段，圓管周圍的流場進行了試驗研究，分析了圓管周圍時均水平流速、時均垂直流速以及渦脫落過程。Smith等[4]分別采用了結合k-ε湍流模型進行封閉的雷諾平均Navier-Stokes(Reynolds average Navier-Stokes,簡稱RANS)模型及大渦模擬(large eddy simulation，簡稱LES)模型對這5個特征階段圓管周圍流場進行了數值模擬，通過與Jensen等[3]的試驗結果進行比較，驗證了兩種模型的可行性。計算結果表明，LES模型對于水平流速及垂直流速的預測更為精準，而結合k-ε模型封閉的RANS方法對于渦尺度的計算準確度受網格尺寸影響較大。

對圓管近水平壁面的研究相對較多。Bearman和Zdravkovich[5]通過試驗探究了圓管海床之間間隙與圓管直徑比值(間隙比)對渦脫落的影響，并分析了圓管周圍的平均壓力系數，試驗表明當該比值等于0.3時，圓管后方出現渦脫落現象，該值也被定義為臨界比值。Zdravkovich[6]通過試驗測量了靠近平坦海床圓管周圍的升力系數及阻力系數，試驗表明升力系數主要由間隙比決定，而阻力系數由邊界層厚度與圓管海床之間比值決定。Lei等[7]的試驗表明了升力系數及阻力系數均由間隙比決定，但受到邊界層影響。Wang和Tan[8]在邊界層厚度與圓管海床之間比值為0.4的情況下進行了試驗，結果表明當間隙比達到0.3時出現了周期性渦脫落。陳鎣等[9]對均勻來流下靠近壁面處在垂直流向做強迫振蕩運動的光滑圓柱水動力特性進行了試驗研究，通過采集順流向和垂直流向的力，得到了阻力系數、升力系數、相位角等與間隙比、振蕩頻率和振幅之間的關系。鄂學全等[10]在試驗中利用平板來模擬海床邊界，測量了圓間隙比對升力系數及阻力系數的影響，試驗結果表明升力系數及阻力系數隨著間隙比增大而減小。周磊等[11]通過粒子圖像測速(particle image velocimetry,簡稱PIV)試驗研究了間隙比對近壁單圓柱繞流的影響，結果表明隨著間隙比的增大，圓柱尾流的流動形態由單個順時針運動的旋渦演變為逐漸對稱的旋渦對，且尺度逐漸減小，同時壁面對圓柱尾流的影響逐漸減弱。在數值模擬方面，Br?rs[12]和Zhao等[13]分別應用k-ω和k-ε湍流模型求解雷諾平均方程，模擬了靠近平坦海床的圓管繞流，獲得的結果與試驗數據相比具有很高的吻合度。謝錦波[14]基于計算流體力學軟件FLUENT，結合k-ω湍流模型對橫管周圍繞流流場進行了模擬計算，計算結果表明所用模型能夠較好地模擬速度場、壓力場和切應力場分布。Ong等[15-17]在高雷諾數情況下應用k-ε模型進行了數值模擬，研究了渦脫落機理。Zhao和Cheng[18]通過數值模擬研究了湍流場中二維近壁圓柱的渦激振動問題，分析了在間隙比分別為0.002和0.3的低間隙比情況下，邊界層對渦激振動的影響。

通過對近壁面圓管繞流問題研究進展的總結可以發現，前人多采用求解RANS方程或者大渦模擬進行計算。考慮到雷諾平均方法在捕捉流場細節方面效果不佳，而大渦模擬又存在對計算資源要求巨大的問題，因此這里采用Spalart-Allmaras改進延遲分離渦模擬(improved detached eddy simulation，簡稱IDDES)以更深入地分析海底管道周圍流場特性，該模型對邊界層內的流動采用雷諾平均方法，對邊界層外的流動采用大渦模擬，在節省計算資源的同時可以模擬湍流的分離，因此可以更深入地分析圓管周圍水動力特征。

1 數學模型

混合模型的誕生主要為了彌補RANS在捕捉流場細節方面的不足，同時避免LES較大的計算量，最初是由Spalart等[19]在1997年提出的一種混合RANS/LES的分離渦(detached eddy simulation，簡稱DES)模型。該模型是一個使用單一湍流模型的三維非定常數值模擬方法，在網格密度滿足大渦模擬要求的區域應用LES模型，在近壁面附近網格密度滿足RANS要求的區域應用標準的Spalart-Allmaras湍流模型[20]解雷諾平均Navier-Stokes方程，從而大幅減少計算時間[21]。

在構造RANS/LES混合模型時，需要對方程中的耗散項進行轉換，方法是構造一個開關函數。RANS/LES的轉換則通過引入長度尺度來判別。DES方法的長度尺度被定義為：

dDES=min(d,ψCDESΔ)

(1)

式中：d為與壁面最近的距離；ψ是低雷諾數修正系數；CDES為經驗常數，建議取值為0.65；Δ為網格尺度。

在壁面附近，d?ψCDESΔ時，采用Spalart-Allmaras湍流模型[19]進行計算；在ψCDESΔ?d區域時，采用亞網格模型進行計算。

但是當網格密度介于LES模型與RANS模型要求之間時，該模型不能捕捉全部的流場紊動，模擬的渦流黏性與雷諾應力將減小，即模擬應力耗盡現象(modeled stress depletion,簡稱MSD)[22]。為了解決由網格密度引起的近壁面處模擬應力耗盡的問題，Spalart等[22]在標準DES模型的基礎上開發了延遲分離渦(delayed detached eddy simulation，簡稱DDES)模型。DDES模型通過渦流黏性檢測邊界層的位置，并強制邊界層內的流場由RANS方法求解。

DDES方法長度尺度為：

dDDES=d-fd·max(0,d-ψCDESΔ)

(2)

式中：fd為與渦流黏性相關的函數。在LES區域(rd?1)，fd等于1，DDES模型即為DES模型；在其余區域fd等于0。

Spalart-Allmaras IDDES模型是由Shur等[23]基于DDES和WMLES(wall modeled LES)改進的RANS/LES混合模型。WMLES模型在邊界層內部區域使用雷諾平均方法，在外部使用大渦模擬，但在RANS和LES模型切換區域存在偏差。IDDES模型解決了該問題，該模型將計算域分為三類子域，分別為遠離壁面區域、近壁面區域、二者之間區域。當入流條件中沒有湍流時，IDDES模型中的DDES方法就會被激活。當入流條件中有湍流且網格密度精細到足夠覆蓋邊界層內的渦時，IDDES模型中的WMLES方法就會被激活。RANS和LES兩種方法通過式(3)結合在一起：

dIDDES=fB(1+fe)dRANS+(1-fB)dLES

(3)

式中：dRANS=d，dLES=ψCDESΔ，fB及fe為經驗函數。

2 計算模型

2.1 計算域及邊界條件設置

利用OpenFOAM模擬了在穩流條件下沖刷初始過程中4個特征階段管道周圍的流動，每一個階段都有固定的沖刷地形。4個沖刷地形取自Mao[2]通過沖刷試驗測得的床型剖面，每個階段的沖坑深度與管道直徑比值G/D分別為0.11、0.42、0.54及0.72。流場計算域總長度選取為40D，計算域高度、寬度分別選取為11D、4D。管道中心放置在(x，y，z)=(0，0.5D，2D)，管道長度為4D，其中D為管道直徑，見圖1。

圖1 模型計算域

圓柱壁面和沖刷床壁面采用無滑移壁面邊界條件，流域上下邊界采用對稱邊界條件，流域兩側邊界采用周期性邊界條件[24]。空間離散采用線性插值格式，時間離散采用向后差分格式。

通過采用雷諾相似準則，設置圓柱直徑D為1 m，入流速度為1 m/s，通過修改黏性系數μ使雷諾數與Jensen等[3]試驗一致(Re=6 000)。

2.2 網格設置

計算網格采用全結構化六面體網格，并對圓管周圍及床面附近的網格進行了加密，如圖2所示。表1列出了每個算例的網格數、近管道壁面及近床面區域的網格厚度，4個算例的流場計算域網格節點數均為1.2×107左右，其中近管道壁面及近床面區域，物面法向方向上第一層計算網格最大厚度均為0.002，y+值均在1左右。

圖2 模型網格

表1 模型網格

2.3 入流邊界層處理

在Jensen等[3]試驗中，邊界層厚度為2D。為了與試驗一致，研究在計算海底管道繞流前，先模擬了穩流在4個無管道的固定床型剖面隨時間充分發展的過程，計算了模型的入流邊界層厚度及速度入口，以驗證邊界層設定與試驗一致。

首先建立了4個無管道床型剖面的計算域及網格，計算域整體長度為40D，高度為11D，寬度為4D。通過OpenFOAM對4個算例進行模擬，至計算域內流動充分發展后達到穩定。計算結果表明，4個模型得到的邊界層厚度與試驗條件十分接近，均為2D左右，表明該方法設定入流邊界條件較為可靠，因此取4個算例穩定后的入流口速度分別作為后續有管道模型數值模擬的入流口速度。

3 模型驗證與結果分析

3.1 模型驗證

Jensen等[3]通過試驗分析了沖刷床剖面上固定管道周圍的流場，試驗中通過聲學多普勒流速儀測量了幾個上游和下游縱向位置的水平和垂直速度剖面。文中研究通過模擬水平管道周圍的流速，與該試驗結果進行比較，以驗證計算模型的可行性及準確性。從模擬結果中提取了與試驗相同監測剖面的速度值，并將試驗數據無量綱化處理后進行對比：

(4)

模型對流速的預測結果通過相對均方根誤差(RMSD)進行客觀評價，均方根誤差定義為：

(5)

圖3和圖4分別為時均水平流速模擬結果及時均垂直流速模擬結果與Jensen試驗數據的對比。圖3表明，4個算例中，數值模擬對于圓柱上游(-3≤x/D≤-1)來流水平流速的預測與試驗結果十分一致，證明了模型定義邊界層方法的可靠性;對圓柱尾流區(0≤x/D≤4)的預測出現了偏差，但偏離程度較低；總體來說，在水平流速的預測上模擬結果較為精準，且可以明顯看出在圓柱后方8D存在擴展邊界層。從圖4可以看出，數值模擬對于垂直流速的預測效果欠佳，4個算例中，對于圓柱上游處(-3≤x/D≤-1)時均垂直速度預測略高于試驗結果，尤其在地形開始變化的區域(x/D=-1)預測偏差逐漸增大，對于圓柱尾流區域(0≤x/D≤4)的計算誤差十分明顯。通過分析，數值模擬對于時均水平流速預測效果很好，但對于時均垂直流速的預測偏差偏大，尤其是在沖刷坑以及沙坡附近的邊界層區域，極有可能是由于沖刷坑深度增加以及沙坡推移，造成邊界層的復雜變化，而IDDES模型在邊界層內區域采用雷諾平均方法求解，再加上復雜邊界的影響使得計算精度降低，過分高估了沙坡的影響，致使對垂直流速的計算高于試驗結果，但整體上流速變化趨勢曲線與試驗曲線接近一致。

圖3 時均水平流速對比(z/D=2)

圖4 時均垂直流速對比(z/D=2)

圖5給出了采用IDDES模型與Smith等[4]采用LES模型和k-ε模型計算得到的時均水平速度均方根誤差對比。在所有算例中，對于圓柱上游(-3≤x/D≤-1)入流邊界層的預測，IDDES模擬結果略顯優勢，平均均方根誤差僅為6.2%，而Smith等[4]在模擬中采用對數方法設定入流邊界層，基于k-ε模型及LES模型4個算例的總平均均方根誤差分別為14.2%及12.3%。對于圓柱尾流區域(0≤x/D≤4)預測，IDDES模型優勢愈發明顯，該模型對于這一區域預測誤差稍有增加，總平均均方根誤差為10.5%。而k-ε模型以及LES模型對于該區域預測誤差大幅度提高，總平均均方根誤差分別為30.0%及21.8%。分析表明采用IDDES模型對于水平流速模擬結果明顯優于LES模型及k-ε模型。

圖5 時均水平速度均方根誤差(z/D=2)

圖6為三者模擬得到的時均垂直速度均方根誤差對比。由于垂直流速較小，均方根誤差值遠低于時均水平速度均方根誤差值。分析表明，在G/D=0.42及G/D=0.54的算例中，對于近尾流區域(1≤x/D≤1.5)，IDDES模型預測誤差雖略有提高但并沒有出現大幅波動，兩個算例總平均均方根誤差為4.2%，而k-ε模型以及LES模型預測誤差出現跳躍式上升，總平均均方根誤差分別為18.7%及16.3%。這可能是由于此時圓柱后方沖刷床出現了坡度較陡的沙坡，使得k-ε模型以及LES模型過度高估了沙坡影響。在其他位置，三者預測結果接近一致。

圖6 時均垂直速度均方根誤差(z/D=2)

Smith等[4]計算采用的兩種方法中，雷諾平均方法由于受到邊界條件、湍流模型條件等影響，其計算結果出現波動性可能較大，對流速預測誤差相對較大；大渦模擬結果優于雷諾方法，但該方法計算三維模型時可靠度最高，計算量也十分龐大，對二維模型的計算仍有欠缺。通過前述分析可以發現IDDES模型彌補了上述兩種方法的不足，計算得到的流速分布與試驗結果吻合度更高，這進一步驗證了該方法的可靠性。

3.2 時均流場分析

將各算例流場的水平流速進行時間平均，匯總如圖7所示。整體來說，尾流場呈現出一定結構性，體現最為明顯的便是圓柱底部的間隙流及圓柱下游回流區，其中間隙流流速受沖刷坑深度影響較大，而回流區范圍與圓柱頂端和底端剪切層及間隙流有關。在G/D=0.11的情況下，圓柱尾流區約10D區域為回流區，該區域中心位置回流流速較大。這種現象的出現可能是由于沖刷坑并未完全發展，深度較淺，從沖刷坑繞過圓柱的水流雖然流速較大，但截面太小，流量仍然有限，大部分流體都是從圓柱上方繞過，此時間隙流強度較弱，因此流體從圓柱上方流過圓柱后產生大范圍回流以補償。而在G/D=0.42的情況下，由于沖刷坑已經發展至一定深度，間隙流強度及范圍變大，圓柱底部出現剪切層，回流區域大幅縮減至1.3D范圍。且由于沙坡的影響，間隙流及回流區均向上偏斜。在G/D=0.54的算例中，隨著沖刷坑發展以及沙坡向后推移，回流區域又逐漸變大至1.7D，間隙流及回流區向上偏移程度減小。在G/D=0.72時，回流范圍擴展至圓柱后方2D處，且回流區域流場逐漸關于圓柱中心所在水平線(y/D=0.5)對稱。

圖7 時均水平流速(z/D=2)

圖8給出了模擬得到的時均垂直流場。可以看到對于圓柱上游，4個算例來流方向幾乎都是水平的，直到流動靠近圓柱時發生分離，形成上升流及下降流。而隨著沖坑的逐漸發展，圓柱前方上升流區域逐漸縮小，而下降流區域逐漸擴大。在圓柱近尾流區均存在上升流，G/D=0.11的情況下由于后方出現的沙坡尺寸較小，上升流范圍很小，但較淺的沖坑致使該區域流速較大。G/D=0.42時，由于沖坑后方沙坡已充分發展至較大尺寸且坡度陡峭，上升流范圍較G/D=0.11的情況明顯增大，二者最大流速相近。在G/D=0.54及G/D=0.72的算例中，隨著沙坡向后推移以及坡度逐漸趨于平緩，上升流區域逐漸縮小，且流速逐漸降低。

圖8 時均垂直流速(z/D=2)

3.3 瞬態流場渦量分析

渦識別技術目前已應用于各個領域，文中采用基于渦量的渦識別方法對圓柱下游渦街情況進行初步分析。在不同工況條件下，取圓柱跨中位置(z=2D)處截面作瞬態渦量圖，如圖9。在G/D=0.11的情況下，由于沖刷坑尚未發展至足夠深度，壁面邊界層對圓柱下游流場影響較大，且從間隙流過的流量有限，間隙流強度較弱。此時圓柱下游近壁面處出現4D長度范圍的低渦量區，且并沒有出現規則的渦脫落現象，僅圓柱頂部出現剪切層脫落，剪切層以小角度向上偏斜后，逐漸出現輕微的波動并向壁面靠近。在G/D=0.42的算例中，間隙流強度增大，圓柱頂部及底部均出現剪切層，底部剪切層在圓柱下游形成了回旋區域，并受到后方沙坡影響向上偏斜與圓柱頂部剪切層接觸融合，有形成渦脫的趨勢。在G/D=0.54的情況下，由于沖刷坑發展且沙坡向后方推移后坡度減緩，壁面邊界層對圓柱底部間隙流影響逐漸減小。此時圓柱下游出現非對稱的雙排尾渦，底端逆時針渦由于壁面邊界層及頂部剪切層的擠壓而受到明顯抑制，尺度小于頂端順時針渦，并最終被抵消至下游僅存單側泄渦。由于沙坡的存在，泄渦向上方偏離，通過沙坡后近乎平行于壁面向下游運動并在圓柱下游8D左右耗散。在G/D=0.72的情況下，由于沖刷坑已經完全展開，壁面邊界層對圓柱底端間隙流影響進一步減小，此時圓柱頂部及底部剪切層在末段逐漸回旋、內扣，圓柱下游出現了穩定且幾乎對稱的雙排渦街，尺度較G/D=0.54條件下的泄渦明顯增大，并大致平行于壁面向下游移動，范圍逐漸擴大。由于此時圓柱后方沙坡已經消失，因此渦脫耗散距離延長至圓柱下游10D左右。

圖9 渦量(z/D=2)

綜合來看，在沖刷前期并無明顯的渦結構產生，中期逐漸產生單側泄渦，最終在后期形成穩定的渦脫落。

4 結 語

基于IDDES方法并利用開源OpenFOAM模擬了沖刷初始過程中4個特征階段在穩流條件下圓管周圍水動力特征，并將模擬結果與Smith等[4]的數值模擬結果以及Jensen等[3]的試驗數據進行比較以驗證模型的適用性，并模擬了管線后方的尾渦形態。主要結論如下：

1)間隙比對近沖刷地形圓柱繞流的流動特性影響顯著。間隙比的不同會影響圓柱與沖刷壁面的相互作用，從而引起尾流特征和旋渦形態的變化。

2)沖刷初期圓柱后方存在大范圍回流區。而隨著間隙比增加，回流區經歷了大幅度縮減又緩慢擴大的過程。

3)在沖刷初期階段，由于受邊壁影響旋渦脫落被抑制；隨著沖刷坑深度逐漸增加，圓柱逐漸脫離邊界層影響，開始產生周期性渦脫落現象，且最終尾流的流動形態由單個旋渦逐漸演變為旋渦對。