低雷諾數下圓柱強迫振動特性光滑粒子流體動力學模擬

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

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

海流繞過海底電纜等管線類結構物流動時，柱體兩側發生交替瀉渦現象，誘發結構物周期性振動，進而導致結構物疲勞損壞。前人研究發現，同等條件下柱體強迫振動和渦激振動時受到的流體力比較接近[1]。因此，開展柱體強迫振動問題研究，對于明確流致渦激振動特性具有重要的借鑒意義[2-3]。

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

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

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

1 SPH數學模型

1.1 控制方程

水體看作是微可壓縮流體，其運動控制方程可表述為下述形式：

(1)

(2)

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

(3)

(4)

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

(5)

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

流體壓力可根據下式由密度變化率轉化得到：

(6)

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

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

(7)

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

(8)

(9)

1.2 邊界條件

如圖1所示，計算域的左右兩側分別設置為入流和出流邊界，通過采用開放邊界算法來實現自由出入流邊界，此算法由Tafuni等[23]提出。在入流邊界和出流邊界處各設置一個粒子緩沖區域，入流和出流緩沖區域內分別布置四層或更多層緩沖粒子。緩沖區粒子和流體粒子互相之間的轉換具體表現為緩沖區的粒子在穿過緩沖區邊界進入流體區域就會變為流體粒子而參與數值計算，反之，流體粒子穿過流體區域邊界進入緩沖區內也會變為緩沖區粒子而不再參與計算。左、右兩側緩沖區內粒子速度采用不同的方法得到，左側緩沖區粒子不設置垂向速度，只給定水平前進方向的速度U，而右側緩沖區粒子的速度則是采用鏡像方式得到，左、右兩側緩沖區內的粒子密度則均是采用鏡像SPH插值獲得。此外，本模型使用了核函數連續性修正算法[24]來解決鏡像區域內局部流體粒子的支持域不完整的問題。

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

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

2 數值分析結果

2.1 計算域設置

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

圖2 計算區域示意Fig. 2 Setup of numerical computational domain

2.2 圓柱繞流

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

表1 升力系數最大值Clmax、斯特勞哈數St和阻力系數均值Cdmean對比

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

圖4 計算域內壓力場分布特征(Re=100)Fig. 4 Pressure field distribution characteristic map in computational domain(Re=100)

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

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

2.3 圓柱的強迫振動

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

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

圖7 圓柱在不同激振頻率下的升阻力系數歷時曲線和升力系數振幅譜(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，雷諾數Re=100時，強迫振動圓柱的鎖定區間如圖8所示，本文計算得到的鎖定區間約為f0/fs=0.75～1.2，Koopmann[6]計算得到的鎖定區間約為f0/fs=0.75～1.25，兩者的計算結果基本趨于一致，差值在誤差允許范圍之內，且二者均在f0/fs=1.0兩側對稱分布。

圖8 A=0.25且Re=100時圓柱強迫振動的鎖定區間Fig. 8 Locking interval of the forced oscillating cylinder at A=0.25， Re=100

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

如圖9所示，在鎖定區間內，以激振頻率比f0/fs=0.9、1.1為例，結構下游的瀉渦模式均呈現出2S模式，在每個振動周期內，圓柱的上下兩側均各有一個渦脫落，這兩個渦的強度一致、方向相反，但是當激振頻率更大時，圓柱后方渦旋方向相同的渦會逐漸融合在一起。

圖9 強迫振動圓柱處于鎖定區間時計算域內的渦量場分布Fig. 9 Vorticity distributions in the computational domain when the oscillating cylinder is in the locked region

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

圖10 圓柱升力系數的最大值Clmax和阻力系數的均值Cdmean隨f0/fs的變化規律(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 結 語

基于GPU并行的DualSPHysic4.4開源代碼，通過使用開邊界算法來實現出入流邊界，建立了基于SPH方法的數值水槽。通過對低雷諾數下圓柱繞流和圓柱強迫振動過程的數值計算，得到了以下結論：

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

2) 流體繞過強迫振動圓柱，在低雷諾數的典型工況下，當雷諾數Re=100時，固定圓柱的振動幅值為A=0.25，強迫振動圓柱的鎖定區間為激振頻率比f0/fs=0.75～1.2，且此時升力系數極值隨激勵頻率比f0/fs的增加而增大，阻力系數均值隨激勵頻率比f0/fs的增加而先增大后減小，圓柱升力和垂向位移之間的相位角隨激勵頻率比f0/fs的增加而減小。