魯小輝， 張 嶔， 楊 博， 梁丙臣,3

(1. 中國(guó)石油化工集團(tuán)有限公司安全監(jiān)管部， 北京 100728；2. 中國(guó)海洋大學(xué)工程學(xué)院， 山東 青島 266100；3. 中國(guó)海洋大學(xué)山東省海洋工程重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266100)

在海洋、海岸工程的實(shí)際應(yīng)用中，海洋結(jié)構(gòu)物中的圓管狀結(jié)構(gòu)物通常成對(duì)出現(xiàn)[1]，如平行布置的海底輸油管線等結(jié)構(gòu)。然而當(dāng)水流經(jīng)過(guò)這些成對(duì)的結(jié)構(gòu)物時(shí)會(huì)因?yàn)榇嬖谖槽E干擾而使流動(dòng)更為復(fù)雜。當(dāng)水流流經(jīng)一對(duì)平行布置的圓管時(shí)，上游圓管產(chǎn)生的尾渦會(huì)和下游的圓管相互作用，并且兩圓管的尾跡受制于兩管的間距影響，從而產(chǎn)生延伸體流態(tài)、重新附著流態(tài)以及共同脫落流態(tài)等復(fù)雜的流態(tài)[2]。當(dāng)兩管之間的間距較近時(shí)，下游圓柱處在上游圓柱的尾渦形成區(qū)域，同時(shí)從上游圓柱表而分離的剪切層會(huì)完全包圍下游圓柱，且剪切層在卷曲產(chǎn)生Karman渦之前不會(huì)附著在下游圓柱表面，此時(shí)尾跡處于延伸體流態(tài)。隨著兩管之間的間距增大，從上游圓柱分離的剪切層會(huì)重新附著到下游的圓柱上，此時(shí)尾跡處于重新附著流態(tài)。而當(dāng)兩管之間的間距進(jìn)一步增大時(shí)，兩圓柱的尾跡都出現(xiàn)Karman渦脫落，此時(shí)尾跡表現(xiàn)為共同渦脫落流態(tài)。上述不同間距條件下的流態(tài)特征同樣還受到入流雷諾數(shù)等因素的影響[3]。

迄今為止，眾多學(xué)者對(duì)于不受壁面影響下的串列雙圓管繞流的特征開展過(guò)實(shí)驗(yàn)研究[3-5]。絕大多數(shù)的實(shí)驗(yàn)研究集中于104量級(jí)的雷諾數(shù)區(qū)間內(nèi)。在數(shù)值模擬方面，Meneghini等[6]利用分?jǐn)?shù)步法模擬了雷諾數(shù)在100～200之間的雙圓管流場(chǎng)相互作用及渦脫落特性。劉松和符松[7]、郭曉玲等[8]、涂佳黃等[9]利用直接數(shù)值模擬方法研究了低雷諾數(shù)下串列雙圓管的流場(chǎng)及渦激振動(dòng)特征。Uzun和Yousuff Hussaini[10], Gopalan和Jaiman[11], 趙偉文和萬(wàn)德成[12]利用RANS和LES混合的方法模擬了雷諾數(shù)為1.66×105條件下不同間距的串列雙圓管流場(chǎng)特征。

對(duì)于近壁面條件下的串列雙圓管繞流的研究目前開展較少，目前多集中于對(duì)水平壁面影響下的流場(chǎng)研究。Harichandan和Roy[13],Bhattacharyya和Dhinakaran[14]開展數(shù)值模擬研究，探究了較低雷諾數(shù)條件下的近壁面串列雙圓管流場(chǎng)特征。Prsic等[15]利用大渦模擬技術(shù)研究了在雷諾數(shù)為1.31×104條件下，兩管間距分別為2和5倍管徑，豎直方向距離水平壁面均為1倍管徑的串列雙圓管繞流場(chǎng)特征。Li[1]在Prsic等[15]的基礎(chǔ)上進(jìn)一步探究了不同豎直方向間距對(duì)流場(chǎng)特征的影響。然而對(duì)于沖刷地形下的串列雙圓管流場(chǎng)，Zhao等[16]基于RANS模型對(duì)不同圓管間距條件下的沖刷過(guò)程開展二維數(shù)值模擬，探究了水平間距對(duì)沖刷深度及水動(dòng)力特征的影響。但是目前對(duì)于沖刷地形下串列雙圓管的三維水動(dòng)力精細(xì)化數(shù)值模擬研究方面較為匱乏。因此本文采用Spalart-Allmaras改進(jìn)延遲分離渦模擬(Improved delayed detached eddy simulation,IDDES)以更深入地研究在不同間距及不同沖刷地形條件下串列雙圓管的流場(chǎng)特征，該模型對(duì)邊界層內(nèi)的流動(dòng)采用雷諾平均方法，對(duì)邊界層外的流動(dòng)采用大渦模擬方法，在節(jié)省計(jì)算資源的同時(shí)可以精確模擬湍流的分離，因此可以更深入地分析圓管周圍流場(chǎng)特性。

1 數(shù)學(xué)模型

1.1 基本方程

對(duì)于不可壓縮黏性流體，在直角坐標(biāo)系下，其運(yùn)動(dòng)規(guī)律可以用Navier-Stokes方程來(lái)描述，連續(xù)性方程和動(dòng)量方程分別為

(1)

(2)

式中：ui為速度矢量；xi為笛卡爾坐標(biāo)系；i=1，2，3；j=1,2,3；ρ為流體密度；p為壓力；t為時(shí)間；v為流體的運(yùn)動(dòng)黏滯系數(shù)。

基于RANS理論，流場(chǎng)中任一點(diǎn)速度及壓力可以分解為：

(3)

(4)

將方程和方程帶入方程和方程中可得：

(5)

(6)

1.2 Spalart-Allmaras DES-IDDES模型

混合模型的誕生主要是為彌補(bǔ)RANS在捕捉流場(chǎng)細(xì)節(jié)方面的不足，同時(shí)混合模型還減少LES較大的計(jì)算量。Spalart等[17]在1997年首次提出了一種混合RANS/LES的分離渦(Detached eddy simulation，DES)模型，該模型通過(guò)單一湍流模型來(lái)實(shí)現(xiàn)對(duì)三維非定常流動(dòng)的精確模擬，在網(wǎng)格密度滿足大渦模擬要求的區(qū)域應(yīng)用LES模型，在壁面附近網(wǎng)格密度滿足RANS要求的區(qū)域應(yīng)用標(biāo)準(zhǔn)的Spalart-Allmaras湍流模型[18]來(lái)解雷諾平均Navior-Storkes方程，從而大幅減少計(jì)算時(shí)間[19]。

在構(gòu)造RANS/LES混合模型時(shí)，需要對(duì)方程中的耗散項(xiàng)進(jìn)行轉(zhuǎn)換，方法是構(gòu)造一個(gè)開關(guān)函數(shù)。RANS/LES的轉(zhuǎn)換則通過(guò)引入長(zhǎng)度尺度來(lái)判別。DES方法的長(zhǎng)度尺度被定義為：

dDES=min(d,ψCDESδ)。

(7)

式中：d為與壁面最近的距離；CDES經(jīng)驗(yàn)常數(shù)，建議取值為0.65；ψ為低雷諾數(shù)校正函數(shù)[20]；δ為網(wǎng)格尺度。

在壁面附近，d?ψCDESδ時(shí)，采用Spalart-Allmaras湍流模型[17]進(jìn)行計(jì)算；與ψCDESδ?d時(shí)，采用亞網(wǎng)格模型進(jìn)行計(jì)算。

但是當(dāng)網(wǎng)格密度介于LES模型與RANS模型要求之間時(shí)，DES模型不能捕捉全部的流場(chǎng)紊動(dòng)，模擬的渦流黏性與雷諾應(yīng)力均將減小，即模擬應(yīng)力耗盡現(xiàn)象(Modeled stress depletion, MSD)[20]。為了解決由網(wǎng)格密度引起的近壁面處模擬應(yīng)力耗盡的問(wèn)題，Spalart等[20]在標(biāo)準(zhǔn)分離渦模擬(Detached eddy simulation,DES)模型的基礎(chǔ)上開發(fā)了延遲分離渦模擬(Delayed detached eddy simulation,DDES)模型。DDES模型通過(guò)渦流黏性檢測(cè)邊界層的位置，并強(qiáng)制邊界層內(nèi)的流場(chǎng)由RANS方法求解。

DDES方法長(zhǎng)度尺度為：

dDDES=d-fdmax(0,d-ψCDESδ)。

(8)

式中fd為基于經(jīng)驗(yàn)公式的遮蔽函數(shù)[20]。在LES區(qū)域(rd?1)時(shí)，fd=1，DDES模型即為DES模型；在其余區(qū)域fd=0。

Spalart-Allmaras IDDES模型是由Shur等[21]基于DDES和WMLES(Wall modeled LES)改進(jìn)的RANS/LES混合模型。WMLES模型在邊界層內(nèi)部區(qū)域使用雷諾平均方法，在外部使用大渦模擬，但在RANS和LES模型切換區(qū)域存在偏差。IDDES模型解決了該問(wèn)題，該模型將計(jì)算域分成三類子域，分別為遠(yuǎn)離壁面區(qū)域、近壁面區(qū)域、二者之間區(qū)域。當(dāng)入流條件中沒(méi)有湍流時(shí)，IDDES模型中的DDES方法就會(huì)被激活。當(dāng)入流條件中有湍流且網(wǎng)格密度精細(xì)到足夠覆蓋邊界層內(nèi)的渦時(shí)，IDDES模型中的WMLES方法就會(huì)被激活。RANS和LES兩種方法通過(guò)以下公式結(jié)合在一起：

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

(9)

式中：dRANS=d；dLES=ψCDESδ；fB和fe為經(jīng)驗(yàn)函數(shù)[21]。

2 計(jì)算模型

2.1 計(jì)算域及邊界條件設(shè)置

本文依據(jù)Zhang[22]在西澳大學(xué)的大型環(huán)形水槽中開展的串列雙圓管渾水沖刷試驗(yàn)，本文依照其試驗(yàn)的流場(chǎng)及沖刷穩(wěn)定后的地形剖面數(shù)據(jù)建立三維數(shù)值水槽，探究不同間距條件下的雙圓管局部流場(chǎng)特征。計(jì)算域如圖1所示，計(jì)算區(qū)域的長(zhǎng)度(x方向)、寬度(z方向)和高度(y方向)分別設(shè)置為60D、4D和11D，其中D為圓管直徑。兩個(gè)平行的圓柱分別放置于對(duì)應(yīng)間距地形的沖刷坑中，G1和G2分別為兩個(gè)圓管正下方?jīng)_刷坑深度，上游圓管1的中心位于距入流邊界(xMin)10D處，下游圓管2中心位于距離入流邊界11D+L處，L為2個(gè)圓管之間的間距。在兩管不同間距的條件下，下游圓管2中心到出流邊界的距離最小為44D，最大為49D。這保證了兩圓管下游尾流的充分發(fā)展。

圖1 計(jì)算域示意圖Fig.1 Schematic of the computational domain

計(jì)算域左側(cè)為入流邊界，其入流流速剖面利用相同的求解程序在一個(gè)未放置圓管的平整底面計(jì)算域中生成。該邊界層入流的生成方法已經(jīng)成功的應(yīng)用于對(duì)沖刷地形下單圓管流場(chǎng)的模擬研究中，并取得了良好的效果(1)曲孟祥, 楊博, 梁丙臣, 等. 基于IDDES方法的沖刷地形海底管道繞流數(shù)值模擬[J]. 海洋工程 (待發(fā)表).。計(jì)算域右側(cè)(xMax)邊界為自由出流邊界，流速ui及壓力p均采用零梯度邊界條件。圓柱壁面和計(jì)算域底部邊界(yMin)為無(wú)滑移壁面邊界，即流速ui=0, 壓力p為零梯度邊界條件。計(jì)算域頂部(yMin)距離圓管足夠遠(yuǎn)，因此在頂部應(yīng)用鋼蓋假定和忽略了對(duì)自由液面變化的模擬。計(jì)算域前(zMin)、后(zMax)邊界均采用周期性邊界條件。計(jì)算域及流場(chǎng)參數(shù)如表1所示。

2.2 地形剖面及網(wǎng)格

圖2給出了不同L/D條件下沖刷地形剖面及計(jì)算網(wǎng)格。由圖2中可以發(fā)現(xiàn)，在L/D≤1.0時(shí),僅在雙圓管下方形成了一個(gè)單一的沖刷坑，其深度隨L/D的增大而增大。當(dāng)L/D>1時(shí)，在兩個(gè)圓管之間形成了一個(gè)較小的沙丘，沙丘的尺度隨著L/D的增大而增大。

表1 入流流場(chǎng)參數(shù)及計(jì)算域幾何特征Table 1 Parameters of the flow field and geometric characteristics of the calculation domain

在L/D=4.9時(shí)，計(jì)算網(wǎng)格采用非結(jié)構(gòu)化六面體網(wǎng)格；在L/D≠4.9時(shí)，計(jì)算網(wǎng)格均采用結(jié)構(gòu)化六面體網(wǎng)格。這是因?yàn)樵贚/D=4.9采用結(jié)構(gòu)化網(wǎng)格較難滿足對(duì)于網(wǎng)格非正交性的要求，從而容易導(dǎo)致計(jì)算發(fā)散。為保證計(jì)算的精度，所有計(jì)算網(wǎng)格均對(duì)圓管周圍及床面附近的網(wǎng)格進(jìn)行了加密(見(jiàn)圖1)。

表2中列出了每個(gè)算例的網(wǎng)格數(shù)、近管道壁面及近床面區(qū)域的網(wǎng)格厚度。所有結(jié)構(gòu)化網(wǎng)格算例的網(wǎng)格總數(shù)在2 000萬(wàn)左右，而非結(jié)構(gòu)化網(wǎng)格總數(shù)略少，約為1 380萬(wàn)。在近管道壁面及近床面區(qū)域，垂直壁面方向上第一層網(wǎng)格最大厚度均為0.002D，以保證y+值均在1左右。這些網(wǎng)格參數(shù)已在先前的研究中證明滿足IDDES模擬要求的經(jīng)度①。

3 計(jì)算結(jié)果與分析

3.1 力系數(shù)分析

圖3給出了在不同L/D的條件下，上、下游圓管的阻力系數(shù)CD和升力系數(shù)CL的變化圖。阻力系數(shù)CD和升力系數(shù)CL定義為：

(10)

(11)

式中：FD和FL分別代表單位長(zhǎng)度圓管所受到的水平和豎直方向的作用力；um為實(shí)驗(yàn)來(lái)流流速。

(橫坐標(biāo)代表無(wú)量綱水平距離，縱坐標(biāo)代表無(wú)量綱豎直距離。The horizontal coordinates represent the dimensionless horizontal distance and the vertical coordinates represent the dimensionless vertical distance)

表2 模型網(wǎng)格參數(shù)Table 2 Mesh parameters

由圖3可以發(fā)現(xiàn)，隨著L/D的增大，上游圓管的阻力系數(shù)略有減小，但其平均值變化不大。隨著L/D的增大，上游圓管的升力系數(shù)逐漸變?yōu)樨?fù)數(shù)，雖然當(dāng)L/D較大時(shí)，瞬時(shí)升力系數(shù)會(huì)大于0，但其平均值仍為負(fù)值。這一點(diǎn)與Li[1]的結(jié)果不同，在水平壁面條件下串列雙圓柱的受力升力系數(shù)始終趨近于0。這表明在沖刷地形下，上游圓管會(huì)受到一個(gè)流場(chǎng)施加的向下的作用力，將上游圓管推向沖刷坑的方向，這更有益于上游管線處于穩(wěn)定狀態(tài)[16]。對(duì)于在L/D較小時(shí)，下游圓管阻力系數(shù)始終為負(fù)值，這表明其受到了流場(chǎng)施加的一個(gè)推力。隨著L/D的增大，下游圓管的阻力系數(shù)逐漸變?yōu)檎怠?/p>

從兩個(gè)圓管的升力系數(shù)的變化過(guò)程可以發(fā)現(xiàn)，在L/D=0.0時(shí)，上、下游圓管的升力系數(shù)的波動(dòng)保持著相同的頻率，這說(shuō)明兩個(gè)圓管的尾跡流態(tài)屬于延伸體流動(dòng)模式，即兩個(gè)圓柱體組合成一個(gè)單一繞流結(jié)構(gòu)[23]。而當(dāng)L/D=0.5或L/D=1.0時(shí)，上游圓管的升力波動(dòng)頻率相對(duì)于下游圓管明顯減弱，這說(shuō)明兩圓管的尾軌模式發(fā)生轉(zhuǎn)換，由上游圓管分離的尾渦受到下游圓管的影響而被抑制。隨著L/D的進(jìn)一步增大，上游圓管的升力波動(dòng)逐漸增加，并且與下游圓管保持近似的頻率，這表明兩圓管的尾軌進(jìn)入共同脫落模式。由于在本研究中入流的雷諾數(shù)較大，尾軌模式發(fā)生轉(zhuǎn)換時(shí)的L/D值會(huì)明顯縮小，這也符合先前的研究成果[3]。

圖3 不同L/D條件下，上、下游圓管力系數(shù)變化Fig.3 Variation of the force coefficients under different L/D conditions

3.2 時(shí)均流場(chǎng)分析

圖4給出了在不同L/D的條件下的時(shí)間平均流線及壓力云圖。從中可以發(fā)現(xiàn)，在L/D=0.0時(shí)，從上游圓管分離的流線完全包絡(luò)了下游的圓管，并且在圓管下游形成了較小的尾軌回流區(qū)。在L/D=0.5時(shí)，可以明顯的看到從上游圓管分離的尾渦連續(xù)的附著在了下游管線的表面，且在兩個(gè)圓管之間形成一個(gè)渦對(duì)[24]。這表明圓管的尾跡流態(tài)進(jìn)入了重新附著模式，也進(jìn)一步驗(yàn)證了圖3中上游圓管的升力發(fā)生變化的原因。此時(shí)受到兩管之間渦對(duì)的影響，從上游圓管分離的流線包絡(luò)范圍略微增大。且兩管之間的下部渦旋迫使流線進(jìn)一步向下偏移，從而使沖刷坑的范圍逐漸增大。當(dāng)L/D=1.0時(shí)，從上游圓管分離的尾渦依然附著在下游管線的表面，但是此時(shí)上部渦旋逐漸消失，而下部渦旋的范圍逐漸增大，形成一個(gè)類似方箱流動(dòng)的流場(chǎng)。類似的現(xiàn)象同樣存在于Li[1]對(duì)水平壁面附近的串列雙圓柱的流場(chǎng)模擬中。該渦旋具有更大的尺度與表面曲率，使得圓管下部的流線進(jìn)一步向下偏移，從而增大沖刷坑的深度，但沖刷坑的水平尺度略有減小。隨著L/D的進(jìn)一步增大兩圓管的尾軌進(jìn)入共同脫落模式，從時(shí)均壓力場(chǎng)變化可以發(fā)現(xiàn)，下游圓管的迎流面的壓強(qiáng)逐漸增大，這說(shuō)明上游圓管的尾流對(duì)下游圓管的影響逐漸減小。從時(shí)均流線變化可以看出，兩個(gè)圓管后部逐漸形成了近乎對(duì)稱的回流區(qū)，且下部的回流區(qū)渦旋相對(duì)于上部更強(qiáng)。由圖2中受力變化可知，上述過(guò)程中伴隨著強(qiáng)烈的尾渦脫落。

3.3 瞬態(tài)渦量場(chǎng)分析

圖5給出了不同L/D條件下瞬時(shí)渦量等值面圖，所選取的渦量的等值面為8.0，其中顏色條代表沿著x方向的瞬時(shí)速度。從圖5中可以發(fā)現(xiàn)，當(dāng)L/D=0時(shí)，在雙圓柱后方形成了一個(gè)單一脫落的尾渦，且形成了較為明顯的肋部渦結(jié)構(gòu)。隨著L/D的增大，尾渦結(jié)構(gòu)逐漸變得復(fù)雜。但是在L/D≤1時(shí)，從圓柱分離的尾渦幾乎不與圓柱下方的床面發(fā)生直接的相互作用，說(shuō)明此時(shí)沖刷坑的發(fā)展主要受制于平均流場(chǎng)的影響。而當(dāng)L/D=2時(shí)，可以清楚的發(fā)現(xiàn)從上游圓管脫落的尾渦沖擊了下游圓管下方的床面，從而增大了下游圓管下方的沖刷坑的深度，導(dǎo)致了在兩管之間形成了較陡的沙丘。隨著L/D的進(jìn)一步增大，從上游圓管脫落的尾渦與下游床面相互作用的位置逐漸變化，從而影響了兩管之間沙丘的位置與形態(tài)。當(dāng)L/D=4.9時(shí)，此時(shí)上游圓管尾渦在未到達(dá)下游圓管之前發(fā)生耗散，從而減弱了下游圓管對(duì)下方?jīng)_刷坑的影響。使得兩管之間沙丘的水平尺度和堆積高度不斷的增高，呈現(xiàn)出與單管沖刷類似的地形特征。

(時(shí)均壓強(qiáng)：Time-averaged pressure. 位于z=2D切面處，橫坐標(biāo)代表無(wú)量綱水平距離，縱坐標(biāo)代表無(wú)量綱豎直距離。Located at z=2D, the horizontal coordinates represent the dimensionless horizontal distance and the vertical coordinates represent the dimensionless vertical distance.)

(等值面為8.0，位于z = 2D切面處。Iso-surface=8.0，Located at z=2D.)圖5 不同L/D條件下瞬時(shí)渦量等值面Fig.5 Iso-surface of instantaneous vorticity for different L/D conditions

4 結(jié)論

(1)隨著L/D的增大，上游圓管的升力系數(shù)逐漸變?yōu)樨?fù)數(shù)。而下游圓管的阻力系數(shù)在L/D較小時(shí)為負(fù)值，并且當(dāng)L/D>1時(shí)突然轉(zhuǎn)為正值。上游圓管的渦脫落頻率在重新附著模式下會(huì)出現(xiàn)明顯的減弱。

(2)在L/D≤1時(shí)，圓管上面的分離渦對(duì)圓管下方床面影響較小，圓管下方的沖刷主要受制于平均水流的影響。且兩管之間渦旋隨著L/D的增大而增大，這進(jìn)一步影響了平均流場(chǎng)與沖刷坑的特征。

(3)在L/D≥2時(shí)，兩圓管的尾跡表現(xiàn)為共同渦脫落流態(tài)，此時(shí)上游圓柱的尾渦將會(huì)對(duì)下游圓柱的沖刷產(chǎn)生較大的影響，從而在兩管之間形成沙丘。