朱敏玲，呂學(xué)強(qiáng)(1. 北京信息科技大學(xué)計(jì)算機(jī)學(xué)院，北京 100101 . 北京信息科技大學(xué) 網(wǎng)絡(luò)文化與數(shù)字傳播北京市重點(diǎn)實(shí)驗(yàn)室，北京 100101)

?

光滑深海立管渦激振動(dòng) DES 模擬

朱敏玲1,2，呂學(xué)強(qiáng)2
(1. 北京信息科技大學(xué)計(jì)算機(jī)學(xué)院，北京 100101 2. 北京信息科技大學(xué) 網(wǎng)絡(luò)文化與數(shù)字傳播北京市重點(diǎn)實(shí)驗(yàn)室，北京 100101)

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

關(guān)鍵詞：DES；深海立管；渦激振動(dòng)；網(wǎng)格尺度；時(shí)間步長(zhǎng)

0　引言

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

1　數(shù)值模型建立

1.1N-S 方程

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

式中，p 為壓強(qiáng)，v為運(yùn)動(dòng)粘性系數(shù)，ρ為密度。

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

渦粘系數(shù)

式中，lk-w為湍流長(zhǎng)度尺度，表達(dá)式為

F1和F2為混合函數(shù)，Pk，Pw為湍流生成項(xiàng)，具體定義根據(jù)參考文獻(xiàn)[7]給出。

1.3DES 方法

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

式中，CDES為網(wǎng)格單元的最大邊長(zhǎng)，CDES=0.65 為常數(shù)。對(duì)通常計(jì)算格式，計(jì)算得到的粘性將過(guò)大，可適當(dāng)減小這個(gè)系數(shù)[8]。

2　研究對(duì)象與數(shù)值方法

2.1計(jì)算域及其離散

以圓的直徑 D 為特征尺度，如圖1 所示，計(jì)算域?yàn)?10 D 25 D。圖中計(jì)算域上游來(lái)流區(qū)域?yàn)?5 D，下游尾流區(qū)域?yàn)?20 D，圓柱距離上、下邊界各為 5 D。文獻(xiàn)[9]表明，該計(jì)算區(qū)域邊界選取對(duì)流場(chǎng)的計(jì)算結(jié)果影響小。

圖1　計(jì)算域示意圖Fig. 1　Fluid computational domain

運(yùn)用 H-O 結(jié)構(gòu)化網(wǎng)格離散計(jì)算域，靠近圓處加密網(wǎng)格，沿徑向逐漸放大，如圖2 所示。

圖2　離散網(wǎng)格及放大圖Fig. 2　Mesh and its local magnification

2.2數(shù)值方法與邊界條件

采用有限體積法（FVM）離散 Navier-Stokes 方程，壓力速度耦合迭代采用 SIMPLE 算法，對(duì)流項(xiàng)采用二階離散格式，擴(kuò)散項(xiàng)采用二階中心差分離散格式。對(duì)動(dòng)量方程、標(biāo)量輸運(yùn)方程采用了欠松弛技術(shù)。計(jì)算精度為 10-5。上下邊界為滑移邊界，圓邊界為無(wú)滑移邊界。流體從左至右流動(dòng)，左側(cè)設(shè)定為速度入口，右側(cè)設(shè)定為壓力出口，壓力參考面為出口面。

3　結(jié)果及分析

3.1網(wǎng)格第一層高度對(duì)結(jié)果影響

文獻(xiàn)[10]給出了 DES 方法在圓柱繞流模擬中的適合網(wǎng)格尺度。隨著 CPU 運(yùn)算能力的提升，可對(duì)網(wǎng)格尺度做更為精細(xì)的研究。文獻(xiàn)[11]按經(jīng)驗(yàn)公式

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

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

表1 中網(wǎng)格數(shù)為 2.4 萬(wàn)，時(shí)間步長(zhǎng) 0.002 s，不同網(wǎng)格第一層高度在 Re = 200 時(shí)的計(jì)算結(jié)果，由表可見，隨著網(wǎng)格第一層高度 Δy 的減小，降低，且前 4 組大于1，在 0.5Δy1之后，小于 1。前 4 組時(shí)均值、幅值隨著Δy 的減小而降低，后 4 組時(shí)均值、幅值隨著 Δy 的減小沒(méi)有明顯降低，而是發(fā)生 2% 以內(nèi)的小幅波動(dòng)。前 4 組 St 為 0.212，后 4 組 St 為 0.216。

表1　不同網(wǎng)格第 1 層高度 Re = 200 計(jì)算結(jié)果Tab. 1　Calculated results of different Δy when Re = 200

圖3 所示為表 2 中計(jì)算結(jié)果隨 Δy 的變化規(guī)律，由圖可見，當(dāng) Δy 減小到使得 y+小于 1 之后，繼續(xù)減小對(duì)時(shí)均值、幅值、St 的影響小于 2%?？梢钥闯?y 處漩渦的典型雷諾數(shù)，也反映粘性的影響隨 y 的變化。DES 方法要求 = 1，0.51得到的最接近 1，且小于 0.51之后對(duì)計(jì)算結(jié)果影響小，因此本文之后的研究中以0.51來(lái)確定第 1 層網(wǎng)格高度。

圖3　Re = 200 計(jì)算結(jié)果隨網(wǎng)格第一層高度變化規(guī)律Fig. 3　Calculated results of different Δy when Re = 200

3.2網(wǎng)格數(shù)量對(duì)結(jié)果影響

圖4　阻力系數(shù)時(shí)均值隨網(wǎng)格數(shù)變化規(guī)律Fig. 4　Drag coefficient of different meshes

圖5　升力系數(shù)幅值隨網(wǎng)格數(shù)變化規(guī)律Fig. 5　lift coefficient of different meshes

圖6　斯特羅哈爾數(shù)隨網(wǎng)格數(shù)變化規(guī)律Fig. 6　Strouhal number of different meshes

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

表2　本文部分網(wǎng)格無(wú)關(guān)解結(jié)果與文獻(xiàn)數(shù)據(jù)的比較Tab. 2　Comparison between calculated results of different meshes with the references

3.3時(shí)間步長(zhǎng)對(duì)結(jié)果影響

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

表3　兩組網(wǎng)格 Re = 200 不同時(shí)間步長(zhǎng)計(jì)算結(jié)果Tab. 3　Calculated results of different time steps when Re = 200

圖7　阻力系數(shù)時(shí)均值隨時(shí)間步長(zhǎng)變化規(guī)律Fig. 7　Drag coefficient of different time steps

圖8　升力系數(shù)幅值隨時(shí)間步長(zhǎng)變化規(guī)律Fig. 8　Lift coefficient of different time steps

圖9　斯特羅哈爾數(shù)隨時(shí)間步長(zhǎng)變化規(guī)律Fig. 9　Strouhal number of different time steps

圖7 為阻力系數(shù)時(shí)均值隨時(shí)間步長(zhǎng)變化規(guī)律，由圖可見，9.7 萬(wàn)網(wǎng)格的阻力系數(shù)時(shí)均值曲線位于 80.0 萬(wàn)網(wǎng)格之下，在相同時(shí)間步長(zhǎng)時(shí)，80.0 萬(wàn)網(wǎng)格的阻力系數(shù)時(shí)均值小于 9.7 網(wǎng)格，更接近文獻(xiàn)中實(shí)驗(yàn)結(jié)果，進(jìn)一步驗(yàn)證了本文 3.2 節(jié)的結(jié)論。升力系數(shù)幅值、斯特羅哈爾數(shù)亦有相同的規(guī)律。

3.4Re = 200 結(jié)果

表4　Re = 200 計(jì)算結(jié)果Tab. 4　Calculated results when Re = 200

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

圖10　Re = 200 瞬時(shí)流場(chǎng)等值線云圖Fig. 10　Transient contours when Re = 200

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

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

4　結(jié)語(yǔ)

本文建立了圓柱繞流渦激振動(dòng)的 CFD 計(jì)算模型，采用 H-O 分塊結(jié)構(gòu)化網(wǎng)格離散計(jì)算域，分析了 DES 方法的網(wǎng)格尺度和時(shí)間步長(zhǎng)，對(duì)圓柱繞流渦激振動(dòng)進(jìn)行了數(shù)值模擬研究。本文研究范圍內(nèi)，可得到如下結(jié)論：

1）基于 SST k-ω 湍流模型的 DES 方法模擬低雷諾數(shù)圓柱繞流渦激振動(dòng)結(jié)果合理；

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

3）增加網(wǎng)格數(shù)量可使計(jì)算結(jié)果更接近實(shí)驗(yàn)結(jié)果，但網(wǎng)格數(shù)量到一定程度后再增加對(duì)結(jié)果改善不明顯。本文中 80.0 萬(wàn)網(wǎng)格結(jié)果最優(yōu)；

4）減小時(shí)間步長(zhǎng)可提高計(jì)算精度，需針對(duì)不同對(duì)象進(jìn)行時(shí)間步長(zhǎng)無(wú)關(guān)解研究。

參考文獻(xiàn)：

[1]胡胡偉. 高聳結(jié)構(gòu)繞流與流固耦合的數(shù)值模擬[D]. 西安: 西安建筑科技大學(xué), 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]曹淑剛, 黃維平, 顧恩凱. 考慮流固耦合的彈性圓柱體渦激振動(dòng)研究[J]. 振動(dòng)與沖擊, 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]蔣仁杰. 圓柱繞流場(chǎng)渦致柱體振動(dòng)的研究[D]. 杭州: 浙江大學(xué), 2013. JIANG Ren-jie. Research on vortex-induced vibrations in the flow around circular cylinders[D]. Hangzhou: Zhejiang University, 2013.

[4]丁林, 張力, 楊仲卿. 高雷諾數(shù)時(shí)分隔板對(duì)圓柱渦致振動(dòng)的影響[J]. 機(jī)械機(jī)械工程學(xué)報(bào), 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-ω 湍流模型的二維圓柱渦激振動(dòng)數(shù)值仿真計(jì)算[J]. 艦船科學(xué)技術(shù), 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ī)葉柵中的應(yīng)用研究[J]. 工程熱物理學(xué)報(bào), 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]時(shí)忠民, 劉名名, 郭曉玲. 計(jì)算域?qū)A柱繞流數(shù)值模擬結(jié)果的影響[J]. 中國(guó)水運(yùn), 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. (未鏈接到本條文獻(xiàn)英文信息)

[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方法數(shù)值模擬圓柱繞流[J]. 艦船科學(xué)技術(shù), 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]魏志理, 孫德軍, 尹協(xié)遠(yuǎn). 圓柱尾跡流場(chǎng)中橫向振蕩翼型繞流的數(shù)值模擬[J]. 水動(dòng)力學(xué)研究與進(jìn)展(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

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

基金項(xiàng)目:北京市教育委員會(huì)科技計(jì)劃科研項(xiàng)目(KM201411232015)、北京市屬高等學(xué)校創(chuàng)新團(tuán)隊(duì)建設(shè)與教師職業(yè)發(fā)展計(jì)劃項(xiàng)目(IDHT20130519)、北京信息科技大學(xué)開放實(shí)驗(yàn)室課題(ICDD2015)資助

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

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

doi：10.3404/j.issn.1672–7619.2016.03.027

中圖分類號(hào)：O353.1；

文獻(xiàn)標(biāo)識(shí)碼：A