光滑深海立管渦激振動 DES 模擬

朱敏玲，呂學強(1. 北京信息科技大學計算機學院，北京 100101 . 北京信息科技大學 網絡文化與數字傳播北京市重點實驗室，北京 100101)

?

光滑深海立管渦激振動 DES 模擬

朱敏玲1,2，呂學強2
(1. 北京信息科技大學計算機學院，北京 100101 2. 北京信息科技大學 網絡文化與數字傳播北京市重點實驗室，北京 100101)

摘要:針對深海立管渦激振動流場建立計算模型，分析了第一層網格高度、網格數量、時間步長對深海立管渦激振動DES模擬的升力系數、阻力系數、斯特羅哈爾數的影響，通過與文獻中實驗、計算數據的對比，說明SST k-ω 湍流模型基礎上的 DES 方法模擬低雷諾數深海立管渦激振動準確合理；網格第 1 層高度對計算精度影響較大，按 0.51確定 DES 方法的第 1 層網格高度可得到滿足要求的。

關鍵詞：DES；深海立管；渦激振動；網格尺度；時間步長

0　引言

在海洋工程中，很多構件為圓柱結構，如鉆井平臺、深海油氣立管等。渦激振動可導致構件疲勞、損壞，影響構件壽命和設備安全。胡衛兵[1]、曹淑剛等[2]以 CFX 為平臺對圓柱繞流渦激振動進行了流固耦合數值模擬，分別對固體結構控制點的位移、壓強和圓柱體的位移響應功率譜、兩向振幅等進行了分析。蔣仁杰[3]運用格子 Bolzmann 方法，研究了彈性支撐單圓柱體渦激振動的力學系統，對圓柱振幅、斯特羅哈爾數與折合速度的關系進行了研究，觀察到了一種主要存在于亞臨界雷諾數區域中的偏移振動形態。丁林等[4]從振動控制的角度研究了分隔板對圓柱渦激振動的影響，圓柱-分隔板結構的振動頻率比單獨的圓柱繞流時低，振動頻率與固有頻率比值保持在 0.4 附近。目前，海洋工程的圓柱繞流渦激振動研究大多針對固體、流固耦合、減振等進行，因此從流場角度研究渦激振動有重要的工程價值。李威[5]等驗證了 SST k-ω 湍流模型對低質量比彈性支撐剛性圓柱體渦激振動問題的有效性。本文采用 SST k-ω 湍流模型[6]基礎上的 DES方法，以海洋工程中深海油氣立管為應用背景模擬渦激振動，分析了第一層網格高度、網格數量、時間步長對立管渦激振動升力系數、阻力系數、斯特羅哈爾數的影響，與文獻結果進行了比較，驗證了數值計算模型的正確性，進一步豐富了海洋工程的理論研究。

1　數值模型建立

1.1N-S 方程

以粘性不可壓空氣為流體介質，不考慮傳熱，N-S方程如下：

式中，p 為壓強，v為運動粘性系數，ρ為密度。

1.2Menter k-ω SST 兩方程湍流模式

渦粘系數

式中，lk-w為湍流長度尺度，表達式為

F1和F2為混合函數，Pk，Pw為湍流生成項，具體定義根據參考文獻[7]給出。

1.3DES 方法

DES 方法用長度尺度替代 k- SST 中長度尺度，從而使得計算區域在附面層使用 k- SST 模型，在主流分離區域使用大渦模擬模型。Δ= max[Δx，Δy，Δz]

式中，CDES為網格單元的最大邊長，CDES=0.65 為常數。對通常計算格式，計算得到的粘性將過大，可適當減小這個系數[8]。

2　研究對象與數值方法

2.1計算域及其離散

以圓的直徑 D 為特征尺度，如圖1 所示，計算域為 10 D 25 D。圖中計算域上游來流區域為 5 D，下游尾流區域為 20 D，圓柱距離上、下邊界各為 5 D。文獻[9]表明，該計算區域邊界選取對流場的計算結果影響小。

圖1　計算域示意圖Fig. 1　Fluid computational domain

運用 H-O 結構化網格離散計算域，靠近圓處加密網格，沿徑向逐漸放大，如圖2 所示。

圖2　離散網格及放大圖Fig. 2　Mesh and its local magnification

2.2數值方法與邊界條件

采用有限體積法（FVM）離散 Navier-Stokes 方程，壓力速度耦合迭代采用 SIMPLE 算法，對流項采用二階離散格式，擴散項采用二階中心差分離散格式。對動量方程、標量輸運方程采用了欠松弛技術。計算精度為 10-5。上下邊界為滑移邊界，圓邊界為無滑移邊界。流體從左至右流動，左側設定為速度入口，右側設定為壓力出口，壓力參考面為出口面。

3　結果及分析

3.1網格第一層高度對結果影響

文獻[10]給出了 DES 方法在圓柱繞流模擬中的適合網格尺度。隨著 CPU 運算能力的提升，可對網格尺度做更為精細的研究。文獻[11]按經驗公式

估算 y+，并確定第一層網格控制點離壁面的距離。本文中，按式（8）估算得到的第一層網格高度稱為 Δy1。

NASA 粘性網格厚度計算器[12]也可估算網格第一層高度。NASA 粘性網格厚度計算器是基于空氣介質在平板湍流中按照 Sutherland’s law[13]來計算空氣粘性厚度，估算Δy。本文中，按 NASA 粘性網格厚度計算器估算得到的第一層網格高度稱為 Δy2。計算得到 Δy1> Δy2。

表1 中網格數為 2.4 萬，時間步長 0.002 s，不同網格第一層高度在 Re = 200 時的計算結果，由表可見，隨著網格第一層高度 Δy 的減小，降低，且前 4 組大于1，在 0.5Δy1之后，小于 1。前 4 組時均值、幅值隨著Δy 的減小而降低，后 4 組時均值、幅值隨著 Δy 的減小沒有明顯降低，而是發生 2% 以內的小幅波動。前 4 組 St 為 0.212，后 4 組 St 為 0.216。

表1　不同網格第 1 層高度 Re = 200 計算結果Tab. 1　Calculated results of different Δy when Re = 200

圖3 所示為表 2 中計算結果隨 Δy 的變化規律，由圖可見，當 Δy 減小到使得 y+小于 1 之后，繼續減小對時均值、幅值、St 的影響小于 2%。可以看成 y 處漩渦的典型雷諾數，也反映粘性的影響隨 y 的變化。DES 方法要求 = 1，0.51得到的最接近 1，且小于 0.51之后對計算結果影響小，因此本文之后的研究中以0.51來確定第 1 層網格高度。

圖3　Re = 200 計算結果隨網格第一層高度變化規律Fig. 3　Calculated results of different Δy when Re = 200

3.2網格數量對結果影響

圖4　阻力系數時均值隨網格數變化規律Fig. 4　Drag coefficient of different meshes

圖5　升力系數幅值隨網格數變化規律Fig. 5　lift coefficient of different meshes

圖6　斯特羅哈爾數隨網格數變化規律Fig. 6　Strouhal number of different meshes

針對 Re = 200 的圓柱繞流流場，設置了網格數從2.4 萬逐漸增加至 160.7 萬的 7 組網格，時間步長均為0.001 s，計算結果隨網格數的變化規律如圖4 ～ 6 所示。隨著網格數量的增加，阻力系數時均值、升力系數幅值、斯特羅哈數均有下降的趨勢。網格數為 2.4萬時，阻力系數時均值為 1.553、升力系數幅值為0.808、斯特羅哈數為 0.216，網格數量增加至 160.7 萬時，阻力系數時均值為 1.398、升力系數幅值為0.517、斯特羅哈數為 0.194。表 2 為本文部分網格無關解結果與文獻數據的比較，由表可知，增加網格數量使得計算結果更接近文獻中實驗結果。80.0 萬網格幅值與文獻[14]相差為 8.4%，St 相差為 8.9%，。160.7 萬網格幅值與文獻[14]相差為 26.1%，St 相差為 2.1%。渦激振動中最主要的性能參數為升力系數和斯特羅哈爾數，綜合考慮 2 個參數的結果，80.0 萬網格與實驗結果更接近。與實驗結果相差的可能原因在于計算模型是二維剛體，實驗中為三維氣動彈性模型。因此，下文將使用 80.0 萬網格數量進行研究。

表2　本文部分網格無關解結果與文獻數據的比較Tab. 2　Comparison between calculated results of different meshes with the references

3.3時間步長對結果影響

針對網格數量為 9.7 萬和 80.0 萬的 2 個網格，分別采取 0.1 s 到 1E-5 s 五組時間步長對 Re = 200 圓柱繞流流場進行計算，得到的結果如表 3、圖7 ～ 圖9 所示，由圖可見，2 組網格的阻力系數時均值、升力系數幅值、斯特羅哈爾數都隨著時間步長的減小而先增大，并在 0.001 s 之后趨于穩定。可知 0.001 s 的時間步長在該網格數量時較為適用，繼續降低時間步長對計算精度提升不大，下文將采用 0.001 s 的時間步長。

表3　兩組網格 Re = 200 不同時間步長計算結果Tab. 3　Calculated results of different time steps when Re = 200

圖7　阻力系數時均值隨時間步長變化規律Fig. 7　Drag coefficient of different time steps

圖8　升力系數幅值隨時間步長變化規律Fig. 8　Lift coefficient of different time steps

圖9　斯特羅哈爾數隨時間步長變化規律Fig. 9　Strouhal number of different time steps

圖7 為阻力系數時均值隨時間步長變化規律，由圖可見，9.7 萬網格的阻力系數時均值曲線位于 80.0 萬網格之下，在相同時間步長時，80.0 萬網格的阻力系數時均值小于 9.7 網格，更接近文獻中實驗結果，進一步驗證了本文 3.2 節的結論。升力系數幅值、斯特羅哈爾數亦有相同的規律。

3.4Re = 200 結果

表4　Re = 200 計算結果Tab. 4　Calculated results when Re = 200

根據本文 3.1 ～ 3.3 節的研究，得到 Re = 200 時圓柱繞流渦激振動 DES 模擬結果，如表 4、圖10 ～ 圖12所示。圖10中可見旋渦交替從圓柱兩側脫落，導致圓柱受到來自流場的阻力、升力均隨時間周期性的脈動，升力脈動頻率由旋渦脫落頻率決定，由表 4 可見，斯特羅哈爾數為 0.207，旋渦脫落頻率 f 為 6.06 Hz。由圖11 和 12 可見，升力脈動頻率為阻力的一半，升力的時均值趨于 0。

圖10　Re = 200 瞬時流場等值線云圖Fig. 10　Transient contours when Re = 200

圖11　Re = 200 阻力時程曲線Fig. 11　Drag coefficient vs time when Re = 200

圖12　Re = 200 升力時程曲線Fig. 12　Lift coefficient vs time when Re = 200

4　結語

本文建立了圓柱繞流渦激振動的 CFD 計算模型，采用 H-O 分塊結構化網格離散計算域，分析了 DES 方法的網格尺度和時間步長，對圓柱繞流渦激振動進行了數值模擬研究。本文研究范圍內，可得到如下結論：

1）基于 SST k-ω 湍流模型的 DES 方法模擬低雷諾數圓柱繞流渦激振動結果合理；

2）隨著網格第一層高度 Δy 的減小，y+降低。按0.5Δy1確定 DES 方法的第一層網格高度可得到滿足要求的；

3）增加網格數量可使計算結果更接近實驗結果，但網格數量到一定程度后再增加對結果改善不明顯。本文中 80.0 萬網格結果最優；

4）減小時間步長可提高計算精度，需針對不同對象進行時間步長無關解研究。

參考文獻：

[1]胡胡偉. 高聳結構繞流與流固耦合的數值模擬[D]. 西安: 西安建筑科技大學, 2010. HU Wei. Numerical simulation of wind pass high-rise structure and fluid-solid coupling[D]. Xi’an: Xi’an University of Architecture and Technology, 2010.

[2]曹淑剛, 黃維平, 顧恩凱. 考慮流固耦合的彈性圓柱體渦激振動研究[J]. 振動與沖擊, 2015, 34(1): 58–62. CAO Shu-gang, HUANG Wei-ping, GU En-kai. Study on vortex-induced vibration of an elastic cylinder considering fluidstructure interaction[J]. Journal of vibration and shock, 2015, 34(1): 58–62.

[3]蔣仁杰. 圓柱繞流場渦致柱體振動的研究[D]. 杭州: 浙江大學, 2013. JIANG Ren-jie. Research on vortex-induced vibrations in the flow around circular cylinders[D]. Hangzhou: Zhejiang University, 2013.

[4]丁林, 張力, 楊仲卿. 高雷諾數時分隔板對圓柱渦致振動的影響[J]. 機械機械工程學報, 2013, 49(14): 133–139. DING Lin, ZHANG Li, YANG Zhong-qing. Effect of splitter plate on vortex-induced vibration of circular cylinder at high Reynolds number[J]. Journal of mechanical engineering, 2013, 49(14): 133–139.

[5]李駿, 李威. 基于SST k-ω 湍流模型的二維圓柱渦激振動數值仿真計算[J]. 艦船科學技術, 2015, 37(2): 30–34. LI Jun, LI Wei. Numerical simulation of vortex-induced vibration of a two-dimensional circular cylinder based on the SST kω turbulent model[J]. Ship science and technology, 2015, 37(2): 30–34.

[6]STRELETS M. Detached eddy simulation of massively separated flows[R]. AIAA-01–0879. San Antonio, Texas: American Institute of Aeronautics &Astronautics, 2001.

[7]MENTER R. Zonal Two equation k- turbulence models for aerodynamics flows[R]. AIAA-93–2906. Orlando, FL: American Institute of Aeronautics &Astronautics, 1993.

[8]顧春偉, 陳美蘭, 李雪松, 等. DES模型在壓氣機葉柵中的應用研究[J]. 工程熱物理學報, 2008, 29(12): 2007–2010. GU Chun-wei, CHEN Mei-lan, LI Xue-song, et al. Application of DES model in the compressor cascade flow[J]. Journal of engineering thermophysics, 2008, 29(12): 2007–2010.

[9]時忠民, 劉名名, 郭曉玲. 計算域對圓柱繞流數值模擬結果的影響[J]. 中國水運, 2013, 13(7): 83–88. SHI Zhong-min, LIU Ming-ming, GUO Xiao-ling. Effect of flow field on flow around circular cylinder[J]. China water transport, 2013, 13(7): 83–88. (未鏈接到本條文獻英文信息)

[10]SPALART P R. Young-person’s guide to detached-eddy simulation grids[R]. NASA/CR-2001–211032. Washington: Boeing Commercial Airplanes, 2001

[11]常書平, 王永生, 龐之洋. 用基于SST模型的DES方法數值模擬圓柱繞流[J]. 艦船科學技術, 2009, 31(2): 30–33. CHANG Shu-ping, WANG Yong-sheng, PANG Zhi-yang. Numerical simulation of flow around circular cylinder using SST DES model[J]. Ship science and technology, 2009, 31(2): 30–33.

[12]NASA. Viscous grid spacing calculator[EB/OL]. 1997. http://geolab.larc.nasa.gov/APPS/YPlus/.

[13]SUTHERLAND W. The viscosity of gases and molecular force[J]. Philosophical magazine, 1893, 36(223): 507–531

[14]BRAZA M, CHASSAING P, MINH H H. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder[J]. Journal of fluid mechanics, 1986, 165: 79–130

[15]魏志理, 孫德軍, 尹協遠. 圓柱尾跡流場中橫向振蕩翼型繞流的數值模擬[J]. 水動力學研究與進展(A輯), 2006, 21(3): 299–308. WEI Zhi-li, SUN De-jun, YIN Xie-yuan. A numerical simulation of flow around a transversely oscillating hydrofoil in the wake of a circular cylinder[J]. Journal of Hydrodynamics, 2006, 21(3): 299–308.

Numerical Simulation of Vortex Induced Vibration on Smooth Marine Riser Using DES Model

ZHU Min-ling, LV Xue-qiang
(1. School of Computer Science, Beijing Information Science and Technology University, Beijing 100101, China 2. Beijing Key Laboratory of Internet Culture and Digital Dissemination Research, Beijing Information Science and Technology University, Beijing 100101, China)

Abstract:DES (Detached Eddy Simulation) method was used to simulate vortex induced vibration of smooth marine riser. The height of first layer of the grid and Grid-independent solution and time step-independence solution is obtained. The lift coefficient, the drag coefficient, Strouhal number (St) and other results agree well with experimental data and other numerical results. The results show that, DES method based on SST k-ω turbulence model is credible and valid to simulate vortex induced vibration of smooth marine riser; the requirement of the first layer of the grid can be satisfied by 0.5Δy1.

Key words:DES；Marine riser；Vortex induced vibration；Grid；Time step

作者簡介:朱敏玲(1979–)，女，博士，講師，主要從事計算機仿真及嵌入式領域研究。

基金項目:北京市教育委員會科技計劃科研項目(KM201411232015)、北京市屬高等學校創新團隊建設與教師職業發展計劃項目(IDHT20130519)、北京信息科技大學開放實驗室課題(ICDD2015)資助

收稿日期:2015–03–30； 修回日期: 2015–07–09

文章編號：1672–7619(2016)03–0128–04

doi：10.3404/j.issn.1672–7619.2016.03.027

中圖分類號：O353.1；

文獻標識碼：A