何建勇，高洋洋，王立忠，沃恩海，張卓先

(1. 浙江大學 海洋學院，浙江 舟山 316021； 2. 浙江省海港投資運營集團有限公司，浙江 寧波 315040)

懸鏈線式細長結構在深海油氣懸鏈線立管、錨泊系統(tǒng)及海上風電管纜等海洋工程中廣泛應用。超細長懸鏈線管體結構由于自重和環(huán)境荷載作用，展向曲率沿程變化，相比于豎直管體，循環(huán)往復的海流誘發(fā)變曲率管體產(chǎn)生渦激振動的過程更為復雜。上部平臺大幅運動下會引起懸鏈線管體結構觸底點位置不斷變化，導致管體傾斜角度和曲率劇烈變化。在渦激振動與平臺運動作用下極易引起懸鏈線管體產(chǎn)生疲勞破壞，導致立管油氣泄漏、錨鏈走錨等事故發(fā)生。

目前關于管體繞流流場特性研究主要針對豎直圓柱和傾斜圓柱，開展了大量試驗研究和數(shù)值模擬。研究表明圓柱流場特性主要和雷諾數(shù)密切相關，隨著雷諾數(shù)的增大，圓柱尾流會從層流轉變?yōu)橥牧?，轉變過程中伴隨著尾流過渡、剪切層過渡和邊界層過渡，最明顯的特征是不穩(wěn)定性。Williamson[1-2]研究發(fā)現(xiàn)圓柱尾流旋渦形成區(qū)域的不穩(wěn)定性導致了尾流的紊亂。該不穩(wěn)定性包括小規(guī)模流向渦的形成，大規(guī)模的旋渦變形以及剪切層的不穩(wěn)定。隨著雷諾數(shù)的增大，旋渦脫落會呈現(xiàn)不連續(xù)性、斑狀大尺度結構等特征，后駐點壓力對尾流的不穩(wěn)定性特征更為敏感[3]。Karniadakis和Triantafyllou[4]研究發(fā)現(xiàn)在雷諾數(shù)Re=200附近尾流特征由二維過渡為三維，隨著雷諾數(shù)增大，靠近圓柱展向渦呈波浪狀，下游渦結構則呈肋條狀。Mittal和Balachandar[5]研究表明二次渦的存在對圓柱表面壓力分布沒有顯著的影響。Behara和Mittal[6]研究在150≤Re≤350區(qū)間內(nèi)圓柱尾流的轉捩特性時發(fā)現(xiàn)，波浪狀的展向渦充分發(fā)展必然引起旋渦錯位，導致渦脫頻率隨時間變化，升阻力系數(shù)均值相應減小。Shlrakashi等[7]試驗研究了傾斜角度變化對圓柱尾流中卡門渦街的干擾作用，當傾角增大時，圓柱軸向出現(xiàn)的二次流會導致渦脫頻率減小，旋渦的規(guī)則性降低。Matsuzaki等[8]通過風洞試驗發(fā)現(xiàn)減小長徑比或增大傾斜角度，圓柱底端脫落的旋渦對整體展向渦的干擾作用增強，沿圓柱軸向渦脫強度增大。Razali等[9]試驗研究了傾斜角度對圓柱尾流速度和渦量三維特征的影響，當傾斜角度為45°時，圓柱展向出現(xiàn)軸向流，尾流中存在大量流向渦，軸向流干擾卡門渦街的形成，導致渦脫變形，尾流的三維特性增強。Zhou等[10]研究發(fā)現(xiàn)當圓柱的傾斜角度大于15°時，尾流中出現(xiàn)二次軸向渦，當傾斜角度增大到45°時，速度頻譜圖中峰值減小，譜峰值區(qū)域變寬。Zhao等[11]通過三維數(shù)值模擬發(fā)現(xiàn)隨著傾斜角度的增大，展向渦形狀和圓柱軸線保持平行，圓周時均壓力系數(shù)和無量綱旋渦脫落頻率減小。Wang等[12]研究發(fā)現(xiàn)隨著傾斜角度增大，展向渦的一致性明顯降低而流向渦一致性增加，在向下游發(fā)展的過程中展向渦之間的距離增大，導致渦脫之間相互作用減弱。

對于懸鏈線管體，沿展向曲率變化較大，其尾流特征比豎直和傾斜管體更為復雜。Miliou等[13]數(shù)值研究了雷諾數(shù)Re為100和500條件下彎管尾流中渦形成及脫落特征，結果表明，彎管凹形布置時，由于彎管展向強烈的軸向流存在，導致尾流中渦脫被完全抑制。Gallardo等[14]采用直接數(shù)值模擬方法研究了Re=3 900下彎管凸形布置時曲率對彎管尾流的影響，結果表明彎管豎直段的渦管形態(tài)與彎管軸線保持平行，而靠近水平段由于局部傾角較大，導致渦脫強度減弱，尾流的三維特性增強。Assi等[15]通過粒子圖像測速技術研究了Re=1 000條件下彎曲管體的尾流形態(tài)，研究表明凹形布置時水平段的存在會使來流出現(xiàn)分離，同時產(chǎn)生的旋渦沖擊彎曲段，造成管體整體尾流湍動性增強，而凸形彎曲管體展向渦脫相關性強，整體渦脫更加規(guī)律。

綜上所述，已有研究主要集中于雷諾數(shù)變化對豎直和傾斜圓柱流場形態(tài)及受力特性的影響研究，對于懸鏈線管體彎曲段尾流流場特征及受力特性研究較少。由于局部傾角沿彎管展向變化，產(chǎn)生的軸向流不僅使尾流中展向渦發(fā)生傾斜，同時會對渦脫起到抑制作用，引起時均壓力系數(shù)和回流區(qū)的變化。采用計算流體力學開源代碼OpenFOAM開展了不同雷諾數(shù)條件下懸鏈線彎曲管體繞流三維流場特性研究，主要分析了雷諾數(shù)對凸形布置彎管瞬時、時均流場特征，表面壓力分布和渦脫頻率變化特性的影響規(guī)律。

1 數(shù)學方法

1.1 計算模型

對于不可壓縮黏性流體，在笛卡爾坐標系中，三維Navier-Stokes(N-S)方程和連續(xù)性方程為：

(1)

(2)

式中:xi為x,y,z三個方向上的坐標；ui為xi方向的速度；t，ρ，p和ν分別表示時間、流體密度、壓力和流體運動黏度。

大渦模擬采用濾波函數(shù)將流場中的旋渦分解為大小尺度兩部分。大尺度結構可以通過直接求解N-S方程得到，小尺度結構用亞格子模型求解。經(jīng)過濾波處理的N-S方程和連續(xù)性方程為：

(3)

(4)

采用有限體積法進行離散。在100≤Re≤500和1 000≤Re≤3 900范圍內(nèi)分別采用直接數(shù)值模擬和大渦模擬方法開展彎管繞流流場三維數(shù)值模擬。在低雷諾數(shù)100≤Re≤500范圍內(nèi)，時間項采用二階Crank-Nicolson隱式離散，對流項、壓力項和拉普拉斯項均采用二階高斯線性格式離散。在高雷諾數(shù)1 000≤Re≤3 900區(qū)間，時間項離散采用backward二階隱式，對流項采用高斯線性格式離散，壓力項和拉普拉斯項采用高斯線性limited格式離散，且均為二階精度。

1.2 邊界條件

將懸鏈線管體簡化成由長度6.0D的豎直段，半徑為12.5D的四分之一圓環(huán)以及10.0D的水平段構成的彎管模型，其中D=0.01 m為彎管豎直段截面直徑，單彎管繞流的計算域如圖1所示，計算域尺寸為32.5D×20.0D×28.5D(長×寬×高)。入口邊界至彎管豎直段軸線的距離為10.0D，底部邊界至水平段軸線的距離為10.0D，兩側邊界距離彎管豎直段軸線為10.0D，彎管的曲率半徑設置為R=12.5D，與Miliou等[13]、Gallardo等[14，16]曲率保持一致，豎直段和水平段的長度分別為6.0D和10.0D，以消除頂部和出口邊界對彎曲段的影響[17-18]。s為彎曲段的弧長，即s=RΦ，其中Φ為從豎直段末端開始沿彎曲段逆時針方向的角度。定義豎直段與彎曲段交接處s/D=0，水平段與彎曲段交接處s/D=19.6。

圖1 單彎管繞流計算域示意Fig. 1 Computational domain of flow past single curved cylinder

入口邊界采用狄利克雷邊界條件，流向速度u分別設置為0.01 m/s、0.03 m/s、0.05 m/s、0.10 m/s、0.15 m/s 和0.39 m/s，對應雷諾數(shù)Re=100、300、500、1 000、1 500和3 900，橫向速度與展向速度設置為v=w=0。出口邊界采用紐曼邊界條件，流速沿出口方向梯度和壓力均設置為0。彎管表面采用無滑移邊界條件。兩側邊界和上下邊界設置為自由滑移邊界條件。當雷諾數(shù)處于100≤Re≤500和1 000≤Re≤3 900時，無量綱時間步長ΔtU/D分別取0.01和0.001，計算過程中最大庫朗數(shù)CFL維持在0.55左右，以保證數(shù)值模擬的穩(wěn)定性和收斂性。

1.3 網(wǎng)格敏感性驗證

首先在Re=500和3 900條件下，采用三套不同密度的網(wǎng)格對豎直圓柱和彎管豎直段開展了網(wǎng)格敏感性研究，以確保后續(xù)模擬結果的準確性和可靠性。豎直圓柱展向高度與彎管豎直段保持一致，即z=6.0D。計算域整體網(wǎng)格和彎管局部網(wǎng)格如圖2所示。

圖2 計算網(wǎng)格示意Fig. 2 Schematic view of the computational mesh

阻力系數(shù)CD、升力系數(shù)CL、壓力系數(shù)Cp和斯特勞哈爾數(shù)St的定義分別為：

(5)

(6)

(7)

(8)

式中:FD和FL分別為圓柱受到的阻力和升力；U為來流速度；A為圓柱沿流向的投影面積，對于豎直圓柱而言，A=LD，L為圓柱展向高度，對于彎管而言，A=(R+D/2+Lv)D，Lv是彎管豎直段的長度；p∞為入口處的壓力大小；fs為旋渦脫落頻率。

表1 Re=500和3 900時網(wǎng)格敏感性研究結果Tab. 1 Results of mesh independence test for Re=500 and 3 900

表2 數(shù)值計算結果與文獻結果對比Tab. 2 Comparisons between the numerical results and previous literature

文中彎管的幾何尺寸與Miliou等[13]一致，將模擬得到的瞬時渦量圖、時均流線圖和時均壓力系數(shù)與文獻[13，29]結果(如圖3～4所示)對比，進一步驗證文中數(shù)值模型的可靠性。如圖3所示，當Re=100時，靠近彎管的渦管軸線垂直，向下游發(fā)展的過程中，渦管發(fā)生彎曲變形，文中模擬結果和Miliou等[13]結果非常吻合。時均流線圖和Canabes[28]的結果基本一致，從圖3(b)中可以發(fā)現(xiàn)Re=100時的回流區(qū)主要集中在彎管軸線后方2D左右的范圍內(nèi)，沿展向回流區(qū)面積逐漸減小。如圖4所示，文中得到的前后駐點線時均壓力系數(shù)變化趨勢和文獻結果一致，即時均壓力系數(shù)絕對值沿展向逐漸減小，該值在彎曲段末端約為-0.1，且Re=500時后駐點線時均壓力系數(shù)絕對值大于Re=100的時均壓力系數(shù)。

圖3 Re=100時彎管尾流特征Fig. 3 Wake flow for a curved cylinder at Re=100

圖4 展向時均壓力系數(shù)Fig. 4 Comparisons of the time averaged pressure coefficients along the spanwise

2 結果與討論

2.1 瞬時流場特征

圖5 不同雷諾數(shù)下彎管的渦量等值面三維圖和俯視圖Fig. 5 Instantaneous iso-surface of vorticity contours of curved cylinder at different Reynolds numbers

圖6 不同雷諾數(shù)下豎直圓柱的渦量等值面三維圖和俯視圖Fig. 6 Instantaneous iso-surface of vorticity contours of straight cylinder at different Reynolds numbers

圖7 Re=100時的瞬時渦量Fig. 7 Instantaneous spanwise vorticity along the span at Re=100

圖8 Re=3 900時的瞬時渦量Fig. 8 Instantaneous spanwise vorticity along the span at Re=3 900

圖9 豎直圓柱瞬時渦量Fig. 9 Instantaneous spanwise vorticity along the span

圖10 不同雷諾數(shù)下不同水平位置處的流向瞬時渦量Fig. 10 Instantaneous streamwise vorticity along the x-direction at different Reynolds numbers

2.2 時均流場特征

圖11為Re=100、500和3 900時沿展向不同位置處的時均流線拓撲，如圖所示，尾流中形成的回流區(qū)關于彎管中軸線對稱，隨著雷諾數(shù)的增大，同一展向位置處的回流區(qū)逐漸減小。當Re=100時，在z/D=27.5，22.5和17.5處的回流區(qū)距彎管圓心2D左右，而Re=500和3 900時，相同位置處的回流區(qū)向下游分別發(fā)展了1.5D和1.0D。由于z/D=12.0位置處局部傾角(Φ=56.44°)較大，回流區(qū)消失。對于彎管而言，沿展向局部傾角增大，靠近彎曲段末端的過程中回流區(qū)逐漸消失。

如圖12所示為不同雷諾數(shù)下彎管中心剖面位置處(y/D=10)的時均流線拓撲。在不同雷諾數(shù)條件下，上游區(qū)域的時均流線變化較為一致，即流線基本平行分布，方向和來流保持一致，流線形狀在接近彎曲段的過程中發(fā)生變化，逐漸與彎管凸面相切。上述變化說明在彎管前駐點線附近產(chǎn)生了軸向流。豎直段后方流線仍然保持水平，但是方向卻與來流相反，表明該范圍內(nèi)出現(xiàn)了回流區(qū)，豎直段附近的回流區(qū)大小隨著雷諾數(shù)的增大逐漸減小。彎曲段流線與來流方向之間存在一定夾角，當100≤Re≤500時，在0°<Φ<42.5°范圍內(nèi)流線與x軸的正方向夾角大于90°，流線存在沿x軸負方向和z軸負方向的分量，說明既有回流現(xiàn)象又有軸向流。當42.5°<Φ<90°時，流線與x軸的正方向夾角小于90°，該區(qū)間只有軸向流，不會出現(xiàn)回流區(qū)。當1 000≤Re≤3 900時，回流區(qū)出現(xiàn)的范圍為0°<Φ<44.5°，相比低雷諾數(shù)有所增大。高洋洋等[32]通過數(shù)值模擬研究了傾斜角度變化對圓柱流場特征的影響，結果表明雷諾數(shù)較小時，若傾斜角度大于15°，圓柱尾流中不會出現(xiàn)回流區(qū)，而在高雷諾數(shù)下，傾斜角度超過30°時回流區(qū)消失。Gallardo等[30]在Re=3 900時觀察到回流區(qū)在Φ≥45°后基本不會出現(xiàn)。這和文中得到的高雷諾數(shù)下彎管尾流中回流區(qū)的分布規(guī)律相似。

圖12 不同雷諾數(shù)下在xoz平面時均流線拓撲Fig. 12 Time-averaged streamline topology at xoz plane for different Reynolds numbers

為了深入探究尾流的展向分布特征，圖13給出了不同z/D處〈u′u′〉/U2和〈v′v′〉/U2沿管體中心剖面位置(y/D=10)的分布。坐標軸中xb定義為22.5-(R-D/2)cosΦ，即對應z在基點線上的橫坐標。從基點到〈u′u′〉/U2最大值位置之間的流向距離即為旋渦形成長度，隨著z/D減小，〈u′u′〉/U2和〈v′v′〉/U2峰值減小，即展向渦脫強度降低，同時向下游移動，表明剪切層的發(fā)展更加穩(wěn)定，旋渦形成長度變大，圖7和8的瞬時展向渦量反映了相同的變化規(guī)律。當Re=100時，在z/D=11.0位置處沒有出現(xiàn)〈u′u′〉/U2峰值，且〈v′v′〉/U2基本為0，由圖7(e)可知管體兩側僅有剪切層存在，無旋渦脫落現(xiàn)象。如圖8(e)所示，當Re=3 900時，在z/D=11.0位置處剪切層卷起產(chǎn)生旋渦脫落，〈u′u′〉/U2出現(xiàn)峰值，〈v′v′〉/U2大于0并出現(xiàn)峰值，但峰值范圍相比于其他位置擴大，能量耗散嚴重，即使旋渦脫落，但渦脫強度仍然較低。隨著雷諾數(shù)的增大，可以發(fā)現(xiàn)相同z/D位置處的〈v′v′〉/U2值增大，即雷諾數(shù)增大，渦脫強度相應增大，和圖5中瞬時三維渦量的變化趨勢相似。

圖13 不同雷諾數(shù)下不同展向位置處的雷諾正應力Fig. 13 The normal Reynolds stress of the streamwise along the centreline of the wake at different depths

2.3 時均壓力系數(shù)

圖14～15為時均壓力系數(shù)沿圓周和展向的分布情況，其中定義前駐點所在的位置處θ=0°，沿圓周順時針方向為正。

圖14 不同雷諾數(shù)下的圓周時均壓力系數(shù)Fig. 14 Time-averaged pressure coefficients along cylinder surface for different Reynolds numbers

如圖14所示，不同展向位置的Cp分布關于θ=180°(后駐點)對稱。在z/D=27.5和22.5位置處的Cp值大小非常接近。在彎曲段上端z/D=17.5處，Cp值相比豎直段有所減小，但變化趨勢仍與豎直段保持一致。在z/D=12.0處局部曲率較大，該處的Cp在θ=0°時明顯小于其他展向位置的值，沿圓周的變化幅值也相應減小。Zhao等[11]對傾斜圓柱時均壓力系數(shù)的研究結果表明，圓周Cp會隨傾斜角度的增大而減小。研究發(fā)現(xiàn)分離點對雷諾數(shù)的變化較為敏感，當100≤Re≤500時，分離點θ≈85°，當1 000≤Re≤3 900時，分離點減小為θ≈75°。

從圖15中可以看到，在不同雷諾數(shù)下，豎直段的前駐點時均壓力系數(shù)Cps基本不會隨展向位置而改變。Cps沿彎曲段逐漸減小，產(chǎn)生的壓力差導致流向發(fā)生變化，因此時均流線在上游區(qū)域靠近彎管時形狀發(fā)生彎曲。相比較前駐點時均壓力系數(shù)，雷諾數(shù)變化對后駐點時均壓力系數(shù)Cpb的影響更為明顯，Cpb的絕對值隨著雷諾數(shù)的增大而增加。另外，彎曲段上部(0°<Φ<45°)的Cpb小于下部(45°<Φ<90°)，因此沿彎管展向回流區(qū)內(nèi)外的壓力差減小，時均流線逆時針方向旋轉，最終與來流方向基本平行。前后駐點時均壓力系數(shù)在彎曲段和水平段交接處(s/D=19.6)二者相等，該位置處時均壓力系數(shù)值約為-0.1，與Miliou等[13]、Gallardo等[30]的研究結果基本一致。

圖15 不同雷諾數(shù)下時均壓力系數(shù)沿展向分布情況Fig. 15 Variation of time-averaged pressure coefficients along the span of the curved cylinder for different Reynolds numbers

2.4 旋渦脫落頻率

為獲得彎管在不同展向位置處的旋渦脫落頻率，對距離彎管軸線后方5D位置處的瞬時橫流向速度進行快速傅里葉變換，開展頻譜分析。圖16表示100≤Re≤3 900的速度頻譜，橫坐標為無量綱的渦脫頻率，縱坐標根據(jù)相應的最大譜峰值進行無因次處理。盡管在不同展向位置處都會出現(xiàn)明顯的頻譜峰值，但是隨著z/D的減小，對應的峰值能量逐漸降低。當Re=300、500和1 000時，頻譜峰值頻率帶隨z/D的減小而變寬，這表明沿展向局部傾角增大導致渦脫強度衰減，影響了脫落旋渦的完整性，能量出現(xiàn)耗散，旋渦的紊亂無序導致譜峰頻率范圍擴大。當Re=1 500和3 900時由于尾流三維特性，頻譜圖中存在較多的諧波。

圖16 不同雷諾數(shù)下不同展向位置處的速度頻譜Fig. 16 Power spectrum of velocity at different depths for different Reynolds numbers

圖17對比了不同雷諾數(shù)條件下彎管在不同展向位置處的無量綱渦脫頻率St，可以看到彎管沿展向的St并非保持一致。當Re=100時，在z/D=27.5、22.5和17.5處的St相等，均為0.167，在z/D=12.0處St略有減小。Miliou[13]在該雷諾數(shù)下觀察到了類似現(xiàn)象，即頻譜能量沿展向衰減，但渦脫頻率保持不變。當300≤Re≤3 900 時，展向St的變化規(guī)律較為相似，即從z/D=12.0到17.5變化逐漸增大，從z/D=17.5到27.5位置處St逐漸減小，在z/D=27.5處的St和豎直管體St值基本保持一致，彎曲段和水平段的存在對豎直段渦脫頻率的影響較小[30-31]。Norberg[33]通過總結豎直圓柱St與Re的關系發(fā)現(xiàn)，當Re≤300時，Re的變化對St的影響較大，而當Re>300時，在亞臨界范圍內(nèi)，Re對St的影響減弱。如圖17所示，當100≤Re≤300時，彎管St的范圍為0.166～0.203，當1 000≤Re≤3 900時，彎管的St基本在0.21左右，說明在高雷諾數(shù)范圍內(nèi)，St對Re的變化并不敏感。由于Re=300正好處在St和Re關系的臨界點，其展向St的變化相對其他雷諾數(shù)變化波動更大。

圖17 不同雷諾數(shù)下不同展向位置處的StFig.17 Variation of St at different depths for different Reynolds numbers

3 結 語

基于計算流體力學開源代碼OpenFOAM開展了不同雷諾數(shù)(100≤Re≤3 900)下，凸形布置彎管繞流三維流場特性研究，分析了彎管瞬時和時均流場特征、時均壓力系數(shù)和旋渦脫落頻率沿彎管展向的變化情況，得到如下結論：

1) 豎直段后方渦管形狀與彎管軸線基本保持平行，渦脫頻率、時均壓力系數(shù)和回流區(qū)范圍不會沿展向變化。由于旋渦脫落的滯后現(xiàn)象，彎曲段渦管開始傾斜，渦脫強度降低，沿展向回流區(qū)范圍減少，時均壓力系數(shù)絕對值降低。

2) 當Re=100時，尾流形態(tài)呈現(xiàn)二維狀態(tài)，水平段不存在剪切層分離和渦脫現(xiàn)象。隨著雷諾數(shù)增加至300≤Re≤500時，流向渦出現(xiàn)，展向渦管呈波浪狀，水平段出現(xiàn)剪切層附著和碎渦。當Re≥1 000時，渦管形狀與彎管軸線基本保持平行，表明展向渦脫同相，水平段附近渦脫無序紊亂，三維特征更為明顯。

3) 隨著雷諾數(shù)的增大，彎管后方回流區(qū)減小。前后駐點時均壓力系數(shù)的絕對值沿彎曲段遞減，導致靠近彎管的流線形狀發(fā)生變化。圓周時均壓力系數(shù)的變化幅度沿展向逐漸減小。相比于前駐點時均壓力系數(shù)，后駐點時均壓力系數(shù)對雷諾數(shù)的變化更為敏感。

4) 低雷諾數(shù)Re=100時彎曲段渦脫頻率和豎直段、水平段基本一致，當300≤Re≤3 900時，彎管展向渦脫頻率存在一定差異，沿彎曲段呈現(xiàn)遞增的趨勢。整體而言，彎管豎直段的渦脫頻率和相同雷諾數(shù)下豎直圓柱的值非常接近。