低雷諾數(shù)下圓柱強(qiáng)迫振動(dòng)特性光滑粒子流體動(dòng)力學(xué)模擬

趙宇蒙，溫鴻杰，任 冰，王 超

(1. 大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2. 中國(guó)建筑股份有限公司，北京 100029；3. 華南理工大學(xué) 土木與交通學(xué)院，廣東 廣州 510641)

海流繞過(guò)海底電纜等管線類結(jié)構(gòu)物流動(dòng)時(shí)，柱體兩側(cè)發(fā)生交替瀉渦現(xiàn)象，誘發(fā)結(jié)構(gòu)物周期性振動(dòng)，進(jìn)而導(dǎo)致結(jié)構(gòu)物疲勞損壞。前人研究發(fā)現(xiàn)，同等條件下柱體強(qiáng)迫振動(dòng)和渦激振動(dòng)時(shí)受到的流體力比較接近[1]。因此，開(kāi)展柱體強(qiáng)迫振動(dòng)問(wèn)題研究，對(duì)于明確流致渦激振動(dòng)特性具有重要的借鑒意義[2-3]。

目前，多位學(xué)者針對(duì)均勻流中的柱體強(qiáng)迫振動(dòng)問(wèn)題開(kāi)展了大量的研究工作。其中，比較典型的工作如：1964年Bishop和Hassan[4]首先試驗(yàn)發(fā)現(xiàn)，當(dāng)外界激勵(lì)頻率近似等于柱體自然瀉渦頻率時(shí)，結(jié)構(gòu)所受到的升力急劇變大，同時(shí)升力與結(jié)構(gòu)位移間的相位角也會(huì)出現(xiàn)躍變現(xiàn)象。Sarpkaya[5]將升力分解為兩個(gè)分量，即慣性力和拖曳力，其中慣性力與振動(dòng)加速度同向，拖曳力與振動(dòng)速度同向。當(dāng)與結(jié)構(gòu)運(yùn)動(dòng)同相位時(shí)，拖曳力為激勵(lì)力，流體向結(jié)構(gòu)輸送能量；反之，則表現(xiàn)為阻尼力，流體消耗結(jié)構(gòu)能量。

在柱體強(qiáng)迫振動(dòng)問(wèn)題的研究中，除了結(jié)構(gòu)受到的升阻力及其與結(jié)構(gòu)位移間的相位差變化特性，振動(dòng)過(guò)程中的頻率鎖定現(xiàn)象以及結(jié)構(gòu)物下游的瀉渦模式，也是研究人員重點(diǎn)關(guān)注的對(duì)象。1967年Koopmann[6]首次試驗(yàn)明確了圓柱體在低雷諾數(shù)情況下的鎖定區(qū)間。Williamson和Roshko[7]進(jìn)一步試驗(yàn)描述了結(jié)構(gòu)后方的瀉渦模式以及尾渦的形成、脫落和耗散全過(guò)程。Placzek等[8]等研究了雷諾數(shù)Re=100時(shí)圓柱強(qiáng)迫振動(dòng)幅值與頻率對(duì)瀉渦模式的影響?；谟邢拊椒?，趙靜[9]研究了雷諾數(shù)Re=100和200以及10 000時(shí)的圓柱強(qiáng)迫振動(dòng)規(guī)律，分析了流體力與結(jié)構(gòu)振動(dòng)頻率間的相位關(guān)系、柱體強(qiáng)迫振動(dòng)的鎖定區(qū)間以及柱體后方的瀉渦模式。梁亮文[10]數(shù)模研究了雷諾數(shù)Re=200時(shí)的圓柱強(qiáng)迫振動(dòng)特性。

由于低雷諾數(shù)下柱體強(qiáng)迫振動(dòng)物模試驗(yàn)對(duì)試驗(yàn)條件要求比較苛刻，實(shí)施難度較大?，F(xiàn)有研究多是采用有限差分、有限體積等基于傳統(tǒng)歐拉網(wǎng)格的數(shù)值模型來(lái)開(kāi)展。這些模型往往需要引入繁瑣的網(wǎng)格重構(gòu)技術(shù)或使用浸入邊界方法，來(lái)近似處理復(fù)雜的物面邊界問(wèn)題。而無(wú)網(wǎng)格光滑粒子流體動(dòng)力學(xué)(SPH)方法，不受網(wǎng)格畸變影響，在處理具有復(fù)雜結(jié)構(gòu)運(yùn)動(dòng)問(wèn)題時(shí)具有很大優(yōu)勢(shì)。因?yàn)橛?jì)算效率和計(jì)算精度的影響，目前采用SPH方法開(kāi)展柱體強(qiáng)迫振動(dòng)的研究文獻(xiàn)不多，如Morris等[11]模擬了低雷諾數(shù)下圓柱繞流的速度場(chǎng)和壓力場(chǎng)；Bouscassea等[12]模擬了雷諾數(shù)為180，佛羅德數(shù)在0.3到2.0之間的近自由表面圓柱繞流問(wèn)題；Sun等[13]及Zhang等[14]研究了鈍體繞流問(wèn)題；Sun等[15]研究了柱體在激勵(lì)頻率為0.8倍其自然渦脫落頻率時(shí)的振動(dòng)特性。近年來(lái)，隨著多種壓力修正算法的提出以及并行計(jì)算技術(shù)的引入，基于SPH方法的水動(dòng)力學(xué)模型的計(jì)算精度和計(jì)算效率均得到了顯著提升。下文借助于DualSPHysic4.4開(kāi)源代碼，研究了激勵(lì)頻率對(duì)圓柱強(qiáng)迫振動(dòng)特性的影響。

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

1.1 控制方程

水體看作是微可壓縮流體，其運(yùn)動(dòng)控制方程可表述為下述形式：

(1)

(2)

式中：ρ，u，P分別是密度、速度和壓力；g是重力加速度；υ0為動(dòng)力黏滯系數(shù)。方程(1)和(2)可寫(xiě)為下述SPH離散形式：

(3)

(4)

Molteni和Colagrossi[16]提出的密度耗散項(xiàng)，可有效減小流體壓力的非物理性震蕩。引入密度耗散項(xiàng)的連續(xù)性方程可以表示為：

(5)

式中：δφ是可變參數(shù)，在文中δφ的值取為0.1[17]，c0為參考聲速。

流體壓力可根據(jù)下式由密度變化率轉(zhuǎn)化得到：

(6)

式中：γ=7，ρ0=1 000 kg/m3。

本模型使用Wendland[18]提出的五次函數(shù)作為核函數(shù)：

(7)

在使用SPH方法模擬強(qiáng)流場(chǎng)問(wèn)題時(shí)，如果流體粒子分布不均，將導(dǎo)致壓力場(chǎng)震蕩，甚至在局部流場(chǎng)出現(xiàn)空腔等非物理性現(xiàn)象。為了增強(qiáng)流場(chǎng)的穩(wěn)定性，Xu等[19]提出了粒子移位矯正法，該方法由Vacondio等[20]引入到了微可壓縮SPH模型中。借鑒Lind等[21]的研究成果，粒子的移位距離δr可以表示為：

(8)

(9)

1.2 邊界條件

如圖1所示，計(jì)算域的左右兩側(cè)分別設(shè)置為入流和出流邊界，通過(guò)采用開(kāi)放邊界算法來(lái)實(shí)現(xiàn)自由出入流邊界，此算法由Tafuni等[23]提出。在入流邊界和出流邊界處各設(shè)置一個(gè)粒子緩沖區(qū)域，入流和出流緩沖區(qū)域內(nèi)分別布置四層或更多層緩沖粒子。緩沖區(qū)粒子和流體粒子互相之間的轉(zhuǎn)換具體表現(xiàn)為緩沖區(qū)的粒子在穿過(guò)緩沖區(qū)邊界進(jìn)入流體區(qū)域就會(huì)變?yōu)榱黧w粒子而參與數(shù)值計(jì)算，反之，流體粒子穿過(guò)流體區(qū)域邊界進(jìn)入緩沖區(qū)內(nèi)也會(huì)變?yōu)榫彌_區(qū)粒子而不再參與計(jì)算。左、右兩側(cè)緩沖區(qū)內(nèi)粒子速度采用不同的方法得到，左側(cè)緩沖區(qū)粒子不設(shè)置垂向速度，只給定水平前進(jìn)方向的速度U，而右側(cè)緩沖區(qū)粒子的速度則是采用鏡像方式得到，左、右兩側(cè)緩沖區(qū)內(nèi)的粒子密度則均是采用鏡像SPH插值獲得。此外，本模型使用了核函數(shù)連續(xù)性修正算法[24]來(lái)解決鏡像區(qū)域內(nèi)局部流體粒子的支持域不完整的問(wèn)題。

圖1 出入流邊界條件示意Fig. 1 Setup of outflow and inflow boundary conditions

上、下邊界采用動(dòng)邊界條件[25]來(lái)模擬，為了減小原動(dòng)力邊界流體粒子與固壁粒子之間的排斥，采用Ren等[26]提出的修正方法對(duì)固壁粒子壓力進(jìn)行修正。本模型中上、下邊界及圓柱分別設(shè)置兩層固避邊界粒子，其物理屬性與流體粒子都是通過(guò)方程(1)和(2)的計(jì)算得到，但是不同情況下的固避邊界粒子的物理屬性又有差別，具體表現(xiàn)在：為了模擬自由滑移邊界條件，上、下邊界中的固避邊界粒子速度給定其水平方向的速度固定為U，不設(shè)置垂向速度，密度值則是通過(guò)將固避粒子鏡像到流體區(qū)域內(nèi)進(jìn)行插值而得到；固定圓柱中固避粒子的位移不更新，其速度值始終為0；強(qiáng)迫振動(dòng)圓柱中的固避粒子速度根據(jù)最初的設(shè)定值保持不變，位移則根據(jù)設(shè)定的速度得到。

2 數(shù)值分析結(jié)果

2.1 計(jì)算域設(shè)置

如圖2所示，本模型中坐標(biāo)原點(diǎn)設(shè)置在圓心處，圓柱的直徑為D=0.1 m，圓柱左邊到入流邊界的距離為10D，計(jì)算域的長(zhǎng)度為60D，寬度為30D，可以足夠保證邊界對(duì)流場(chǎng)無(wú)明顯影響。計(jì)算域的左、右兩側(cè)分別為入流邊界和出流邊界，其緩沖區(qū)內(nèi)均設(shè)置為4層緩沖粒子。本模型在出口邊界前方設(shè)置了垂向的人工阻尼層，寬度為5D，其作用在于通過(guò)減小流體粒子穿過(guò)右側(cè)邊界時(shí)的渦旋強(qiáng)度，來(lái)避免出口處強(qiáng)渦旋對(duì)計(jì)算域內(nèi)壓力場(chǎng)的影響。本模型中初始的粒子間距為Δx=D/50，固避粒子和流體粒子的總數(shù)之和超過(guò)4.6×106個(gè)。在計(jì)算開(kāi)始前，給定流體粒子的初始水平速度為1 m/s，不設(shè)置初始垂向速度，初始密度值為ρ=103kg/m3。雷諾數(shù)的計(jì)算式為Re=UDυ0-1，不同雷諾數(shù)下的固定圓柱和振動(dòng)圓柱的受力特性通過(guò)控制流場(chǎng)速度值不變，圓柱直徑大小不變，只改變黏性值υ0來(lái)實(shí)現(xiàn)。每一工況的模擬時(shí)間均為30 s，運(yùn)算時(shí)間在15.5 h左右，計(jì)算顯卡型號(hào)為NVIDIA GeForce RTX 2080Ti。

圖2 計(jì)算區(qū)域示意Fig. 2 Setup of numerical computational domain

2.2 圓柱繞流

圓柱繞流計(jì)算域如圖2所示。如表1所示，與前人的計(jì)算結(jié)果相比，本文計(jì)算得到的升力系數(shù)最值Clmax和斯特勞哈數(shù)St均吻合較好(St=fsD/U，fs為自然渦脫落頻率，Ts為固定圓柱繞流的瀉渦周期，即Ts=1/fs)，但是阻力系數(shù)平均值Cdmean小于其他人的計(jì)算結(jié)果，Marrone等[27]也發(fā)現(xiàn)在采用與動(dòng)力學(xué)邊界相似的虛擬粒子邊界時(shí)，圓柱的阻力系數(shù)會(huì)偏小。在粒子間距dp分別為D/67、D/50時(shí)，本文的計(jì)算結(jié)果變化非常小，說(shuō)明粒子分辨率已經(jīng)收斂。

表1 升力系數(shù)最大值Clmax、斯特勞哈數(shù)St和阻力系數(shù)均值Cdmean對(duì)比

圖3 升力系數(shù)Cl和阻力系數(shù)Cd歷時(shí)Fig. 3 Time series of the lift and drag coefficients

圖4 計(jì)算域內(nèi)壓力場(chǎng)分布特征(Re=100)Fig. 4 Pressure field distribution characteristic map in computational domain(Re=100)

圖5 右側(cè)出口邊界二個(gè)壓力測(cè)點(diǎn)處壓力歷時(shí)曲線(Re=100)Fig. 5 Time series of the pressure at two measuring points on the outlet boundary(Re=100)

圖6 一個(gè)周期內(nèi)渦的形成、發(fā)展和脫落過(guò)程(Re=200)Fig. 6 The formation, development and shedding of vortices in a period

2.3 圓柱的強(qiáng)迫振動(dòng)

通過(guò)模擬圓柱在雷諾數(shù)Re=100時(shí)的強(qiáng)迫振動(dòng)問(wèn)題，研究圓柱在不同激振頻率下的尾渦脫落模式及其受力特性。文中，f0表示為激振頻率，fv表示為升力頻率，fs表示為自然渦脫落頻率，圓柱振幅為A=ymax/D，其中ymax為圓柱的激振幅值。通過(guò)改變圓柱的激振頻率f0模擬了當(dāng)f0/fs=0.65～1.35之間共計(jì)8種工況下的圓柱強(qiáng)迫振動(dòng)，當(dāng)fv脫離fs，并且鎖定在f0周圍時(shí)，可以認(rèn)為此時(shí)的圓柱位于鎖定區(qū)間之內(nèi)。

不同激振頻率下的升阻力系數(shù)歷時(shí)曲線及升力系數(shù)的振幅譜分別如圖7的左圖和右圖所示。如圖7(a)左圖所示，即當(dāng)f0/fs=0.65時(shí)，升力系數(shù)的拍頻特征明顯；從圖7(a)右圖也可明顯看到兩個(gè)峰值的結(jié)果，此時(shí)圓柱的強(qiáng)迫振動(dòng)位于鎖定區(qū)間之外。升力的主要頻率成分有兩個(gè)：一個(gè)靠近f0，另一個(gè)是fs，此時(shí)圓柱的升力同時(shí)由f0和fs控制，但fs起主要作用。當(dāng)f0/fs=0.75時(shí)，從圖7(b)右圖可以看到升力已經(jīng)主要由fs控制變?yōu)橛蒮0控制，此時(shí)圓柱強(qiáng)迫振動(dòng)已經(jīng)位于鎖定區(qū)間邊界附近，當(dāng)f0/fs=0.8、0.9、1.0、1.1時(shí)，從圖7(c)、(d)、(e)、(f)右圖可以觀察到升力系數(shù)的振幅譜均只有一個(gè)峰值，升力系數(shù)的頻率脫離fs，基本鎖定于f0，圓柱的升力由f0控制，圓柱的強(qiáng)迫振動(dòng)位于鎖定區(qū)間。當(dāng)f0/fs=1.2時(shí)，從圖7(g)左圖可以看到升力中也出現(xiàn)了拍頻現(xiàn)象，右圖的振幅譜也出現(xiàn)了兩個(gè)峰值，但f0仍占主要成分，此時(shí)，強(qiáng)迫振動(dòng)的圓柱即將脫離鎖定區(qū)間。當(dāng)f0/fs=1.35時(shí)，從圖7(h)左圖可以看出升力再次呈現(xiàn)拍頻特征，圖7(h)右圖的振幅譜圖也可明顯觀察到兩個(gè)峰值的結(jié)果，此時(shí)圓柱的強(qiáng)迫振動(dòng)已經(jīng)超出其鎖定范圍。

圖7 圓柱在不同激振頻率下的升阻力系數(shù)歷時(shí)曲線和升力系數(shù)振幅譜(A=0.25， Re=100)Fig. 7 The lift and drag coefficients duration curve and lift coefficient amplitude spectrum of the cylinder under different excitation frequencies at Re=100, A=0.25

在振幅A=0.25，雷諾數(shù)Re=100時(shí)，強(qiáng)迫振動(dòng)圓柱的鎖定區(qū)間如圖8所示，本文計(jì)算得到的鎖定區(qū)間約為f0/fs=0.75～1.2，Koopmann[6]計(jì)算得到的鎖定區(qū)間約為f0/fs=0.75～1.25，兩者的計(jì)算結(jié)果基本趨于一致，差值在誤差允許范圍之內(nèi)，且二者均在f0/fs=1.0兩側(cè)對(duì)稱分布。

圖8 A=0.25且Re=100時(shí)圓柱強(qiáng)迫振動(dòng)的鎖定區(qū)間Fig. 8 Locking interval of the forced oscillating cylinder at A=0.25， Re=100

經(jīng)過(guò)計(jì)算發(fā)現(xiàn)，當(dāng)f0/fs=0.8、0.9、1.0和1.1時(shí)，即在鎖定區(qū)間內(nèi)，圓柱的升力和垂向位移之間的實(shí)時(shí)相位角φ平均值分別為179.2°、123.6°、96.4°、59.2°。Placzek等[6]在當(dāng)f0/fs=0.9、1.0和1.1時(shí)得到的相位角φ分別為113.8°、80.5°、47.7°，經(jīng)過(guò)對(duì)比分析可知，兩者中相位角φ的變化規(guī)律基本一致，即隨著激振頻率的減小而增大。

如圖9所示，在鎖定區(qū)間內(nèi)，以激振頻率比f(wàn)0/fs=0.9、1.1為例，結(jié)構(gòu)下游的瀉渦模式均呈現(xiàn)出2S模式，在每個(gè)振動(dòng)周期內(nèi)，圓柱的上下兩側(cè)均各有一個(gè)渦脫落，這兩個(gè)渦的強(qiáng)度一致、方向相反，但是當(dāng)激振頻率更大時(shí)，圓柱后方渦旋方向相同的渦會(huì)逐漸融合在一起。

圖9 強(qiáng)迫振動(dòng)圓柱處于鎖定區(qū)間時(shí)計(jì)算域內(nèi)的渦量場(chǎng)分布Fig. 9 Vorticity distributions in the computational domain when the oscillating cylinder is in the locked region

如圖10所示，當(dāng)振幅A=0.25且雷諾數(shù)Re=100時(shí)，升力系數(shù)的最大值Clmax隨激振頻率的變化規(guī)律與Placzek等[8]的計(jì)算結(jié)果基本趨于一致，均為先隨著激振頻率的增加先減小而后逐漸增大；阻力系數(shù)均值隨激振頻率的變化規(guī)律也基本相同，二者的均值Cdmean隨激振頻率的增加也均為先增加而后減小，在鎖定區(qū)間左邊界和右邊界都取得極小值后再逐漸增大，雖然變化趨勢(shì)相近，但是本文計(jì)算得到的阻力系數(shù)相比要更小，這是因?yàn)閿?shù)值算法本身造成的差異，不影響對(duì)客觀規(guī)律的揭示。

圖10 圓柱升力系數(shù)的最大值Clmax和阻力系數(shù)的均值Cdmean隨f0/fs的變化規(guī)律(A=0.25， Re=100)Fig. 10 Evolution of Clmax and Cdmean with the frequency ratio f0/fs at A=0.25， Re=100

3 結(jié) 語(yǔ)

基于GPU并行的DualSPHysic4.4開(kāi)源代碼，通過(guò)使用開(kāi)邊界算法來(lái)實(shí)現(xiàn)出入流邊界，建立了基于SPH方法的數(shù)值水槽。通過(guò)對(duì)低雷諾數(shù)下圓柱繞流和圓柱強(qiáng)迫振動(dòng)過(guò)程的數(shù)值計(jì)算，得到了以下結(jié)論：

1) 流體繞過(guò)固定圓柱，在低雷諾數(shù)的典型工況下，當(dāng)雷諾數(shù)Re=100和200時(shí)，通過(guò)數(shù)值模擬得到的升阻力系數(shù)以及斯特勞哈數(shù)，與前人使用其他數(shù)值方法模擬得到的計(jì)算結(jié)果基本一致，顯示SPH數(shù)值模型可以模擬低雷諾數(shù)下流體繞過(guò)固定圓柱后的圓柱受力特性及流場(chǎng)特征。

2) 流體繞過(guò)強(qiáng)迫振動(dòng)圓柱，在低雷諾數(shù)的典型工況下，當(dāng)雷諾數(shù)Re=100時(shí)，固定圓柱的振動(dòng)幅值為A=0.25，強(qiáng)迫振動(dòng)圓柱的鎖定區(qū)間為激振頻率比f(wàn)0/fs=0.75～1.2，且此時(shí)升力系數(shù)極值隨激勵(lì)頻率比f(wàn)0/fs的增加而增大，阻力系數(shù)均值隨激勵(lì)頻率比f(wàn)0/fs的增加而先增大后減小，圓柱升力和垂向位移之間的相位角隨激勵(lì)頻率比f(wàn)0/fs的增加而減小。