曹 興， 劉宇飛， 于 恒， 孔祥鑫， 李慶領(lǐng)

(青島科技大學(xué)機(jī)電工程學(xué)院，青島 266061)

圓柱繞流現(xiàn)象廣泛存在于各類實際工程問題之中[1-2]，如海洋工程、測量設(shè)備、換熱器等。前人對圓柱繞流問題的研究[3-5]主要集中于來流定常流動，但是流體的實際流動過程中存在強(qiáng)烈的非定常性，會產(chǎn)生不同程度的脈動流，如果忽略脈動流對流場的影響，會令研究偏離實際情況。此外，脈動流常用來改變流場與傳熱特性[6]。宋德福等[7]采用數(shù)值方法研究了脈動流下平板壁面的換熱；李國能等[8]采用數(shù)值模擬方法，研究了平行圓柱體在層流脈動流中的相關(guān)特性。但是脈動流條件下的流場中，不僅存在著脈動阻力，而且在垂直來流的展向上會產(chǎn)生較大的升力，同時脈動流在空間上引入了脈動振幅，在時間上引入了脈動頻率，導(dǎo)致脈動流下的流場特性非常復(fù)雜，使得目前少有對于脈動流條件下圓柱繞流特性的研究[9]。相關(guān)研究中，王燕珍[10]對脈動流下圓柱繞流問題進(jìn)行了研究，發(fā)現(xiàn)柱體所受的應(yīng)力在脈動流條件下更大；Armstrong等[11]在湍流脈動條件下研究了圓柱繞流問題，發(fā)現(xiàn)橫向分離的減小與旋渦強(qiáng)度的增大影響了鎖定范圍，但其均未能綜合給出脈動頻率和無量綱脈動振幅對圓柱繞流尾渦脫落、升阻力系數(shù)、旋渦脫落頻率等的影響規(guī)律。為此采用數(shù)值模擬的方法對雷諾數(shù)Re=200、脈動流條件下的圓柱繞流特性進(jìn)行研究，并與相關(guān)文獻(xiàn)進(jìn)行對比，分析流體脈動頻率與振幅對升、阻力系數(shù)以及旋渦脫落頻率的影響規(guī)律，并對尾渦脫落過程的渦量分布進(jìn)行研究。

1 數(shù)值模擬方法

1.1 物理模型與控制方程

計算區(qū)域物理模型如圖1所示。由圖1可知，計算區(qū)域的幾何尺寸為20D×30D(D為圓柱直徑，m，D=0.1 m)，流體入口邊界與圓柱中心的垂直距離為10D，出口邊界與圓柱中心的垂直距離為20D，垂直于流向的上下邊界到圓柱中心的距離各為10D。

圖1 計算區(qū)域物理模型Fig.1 Physical model of computational domain

所采用的控制方程為

(1)

(2)

(3)

式中：μ為流體黏度，(N·s)/m2；ρ為流體密度，kg/m3；u為x方向的速度，m/s；v為y方向的速度，m/s；p為壓力，Pa。

雷諾數(shù)表征流體慣性力和黏性力之比，公式如式(4)所示：

(4)

式(4)中:U為來流速度，m/s。

以水為流動工質(zhì)，數(shù)值計算選用層流模型，Re=200。

1.2 邊界條件與數(shù)值計算方法

對物理模型的邊界條件進(jìn)行了簡化。定義入口為速度入口邊界，給定入口速度呈正弦脈動變化，即：

u(t)=um[1+Asin(2πft)]

(5)

式(5)中：u(t)為入口脈動流速度，m/s；t為流動時間，s;um為穩(wěn)態(tài)流速，m/s；f為脈動頻率，Hz；A為無量綱脈動振幅，A=A0/D，其中A0為脈動振幅，m。

通過改變脈動頻率f和無量綱脈動振幅A(脈動振幅A0與圓柱直徑D的比值，A=A0/D)，來改變?nèi)肟诿}動流工況；定義出口為壓力出口邊界，給定靜壓和適當(dāng)?shù)幕亓鳁l件；圓柱表面設(shè)為無滑移不可滲透固體壁面，近壁面處采用標(biāo)準(zhǔn)壁面函數(shù)法處理；垂直于流向的上下兩邊界采用對稱邊界條件，該邊界上各單元節(jié)點(diǎn)的變量沿法向的分量為零。

對于壓力-速度的耦合選擇Coupled算法，動量方程離散采用Quick格式，對于壓力項離散采用Standard格式，瞬態(tài)項采用二階隱式格式，時間步長設(shè)置為0.002 s。當(dāng)計算域內(nèi)所有控制體積的各方程平均絕對殘差小于10-5時，認(rèn)為迭代計算收斂。

1.3 網(wǎng)格劃分與獨(dú)立性驗證

計算區(qū)域采用結(jié)構(gòu)化網(wǎng)格劃分，出于節(jié)省計算資源同時提高計算精度的考慮，對圓柱周圍區(qū)域的網(wǎng)格采取加密處理，而遠(yuǎn)場區(qū)域的網(wǎng)格相對稀疏。

圓柱水動力可以分解為沿流向的阻力和沿法向的升力，將升力Fl和阻力Fd無量綱化，可以分別得到升力系數(shù)Cl和阻力系數(shù)Cd。

(6)

(7)

式中：Fl為作用于單位長度圓柱上的升力，N/m；Fd為作用于單位長度圓柱上的阻力，N/m。

利用斯特勞哈爾數(shù)，可以表示旋渦脫落的快慢程度。

(8)

式(8)中：fv為旋渦脫落頻率，Hz；U為來流速度，m/s。

當(dāng)Re=200時，選用4套不同密度的網(wǎng)格進(jìn)行獨(dú)立性測試，由4套網(wǎng)格計算得到Cd時均值、Cl振幅、St，如表1所示。由表1可知，隨著網(wǎng)格單元數(shù)的增加，Cd時均值、Cl振幅和St均逐漸增大，且變化的幅度均逐漸變緩。由網(wǎng)格2計算出的Cd時均值、Cl振幅和St與網(wǎng)格1計算出的分別相差1.67%、4.17%和5.26%，網(wǎng)格3的計算結(jié)果與網(wǎng)格2的計算結(jié)果Cd時均值、Cl振幅和St分別相差1.64%、2%和5%，網(wǎng)格4的計算結(jié)果與網(wǎng)格3的計算結(jié)果Cd時均值、Cl振幅和St分別相差0.81%、0.98%和0.95%。當(dāng)網(wǎng)格單元數(shù)大于網(wǎng)格3的單元數(shù)以后，數(shù)值計算結(jié)果已經(jīng)相當(dāng)接近，誤差在1.0%以內(nèi)，計算結(jié)果已經(jīng)滿足網(wǎng)格獨(dú)立性的要求。綜合考慮后選用網(wǎng)格3進(jìn)行計算。網(wǎng)格劃分示意如圖2所示。

表1 網(wǎng)格獨(dú)立性驗證Table 1 Summary of grids independence checks for simulation

圖2 計算區(qū)域網(wǎng)格Fig.2 Schematic of the computational grids

1.4 數(shù)值計算方法驗證

基于文獻(xiàn)[12-14]的計算結(jié)果，對方法的可靠性與準(zhǔn)確性進(jìn)行驗證，結(jié)果如表2所示。在定常流條件下，計算得到的Cd時均值與文獻(xiàn)[12-14]偏差為5.34%～7.46%，Cl振幅與文獻(xiàn)[12-14]偏差為21.54%～31.08%，St的數(shù)值與文獻(xiàn)[12-14]偏差為5%～16.7%， 偏差均在合理范圍之內(nèi)，證明了本文方法的可靠性。

對9種脈動流工況下圓柱繞流場進(jìn)行數(shù)值計算，結(jié)果如表3所示。

表2 計算結(jié)果與文獻(xiàn)結(jié)果對比Table 2 Comparisons between calculated results and literature results

表3 脈動流工況Table 3 The working conditions of pulsating flow

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

2.1 典型時刻渦量分布圖

圖3為不同脈動流工況下，旋渦脫落周期T內(nèi)等時間間隔的4個典型時刻瞬時渦量圖。

圖3 不同脈動流條件下圓柱繞流瞬時渦量Fig.3 Vorticity contours of pulsating flowing around a circular cylinder under different conditions

由圖3可以看出，在t=0時刻，下渦占據(jù)主導(dǎo)地位，分離區(qū)域在圓柱上側(cè)近壁面處形成，此時上個旋渦脫落周期形成渦的剪切層還與近壁面黏連；在t=T/4時刻，下渦占據(jù)主導(dǎo)地位，新的上渦形成，在下渦的拉拽作用下剪切層斷裂；在t=T/2時刻，上渦逐漸開始混合并重新占據(jù)主導(dǎo)地位，下渦加速向下游移動；在t=3T/4時刻，上渦混合成一個大渦并占據(jù)主導(dǎo)地位，在上渦的作用下，新的下渦形成后剪切層被撕扯斷裂離開圓柱表面，此時1個新的旋渦脫落周期完成，如此周而復(fù)始。與定常流條件下圓柱繞流流動分離規(guī)律不同，脈動流條件下圓柱繞流旋渦的剪切層更薄，分離更快，并且由于速度呈脈動變化，脈動流條件下圓柱繞流的渦量圖在1個完整周期內(nèi)沒有明顯的對稱性。

提高脈動頻率f使得尾渦長度變短，旋渦脫落速度加快，混合速度減慢，隨著無量綱脈動振幅A逐漸增大，旋渦生成區(qū)域越來越靠近圓柱表面，這使得處于主導(dǎo)地位渦的脫落速度減緩，另一側(cè)渦的脫落速度加快。

2.2 脈動頻率與脈動振幅對升力系數(shù)的影響

圖4為無量綱脈動振幅A=1，脈動頻率f分別為4、8、10 Hz時升力系數(shù)Cl的時間歷程曲線。對比圖4可以看出，當(dāng)A不變時，隨著f逐漸升高，Cl的振幅逐漸增大。圖5為脈動頻率f=10 Hz，無量綱脈動振幅A分別為0.2、0.6、1時Cl的時間程曲線。對比圖5可以看出，當(dāng)f不變時，Cl的振幅隨著A的增大逐漸增大。與定常流條件下圓柱繞流情況[15]相比可以發(fā)現(xiàn)，在相同雷諾數(shù)下，脈動流條件下Cl的振幅大于定常流條件下Cl的振幅，說明脈動流的作用使圓柱橫向振動增強(qiáng)。

2.3 脈動頻率與脈動振幅對阻力系數(shù)的影響

當(dāng)無量綱脈動振幅A=1時，脈動頻率f分別為4、8、10 Hz時的阻力系數(shù)Cd時間歷程曲線如圖6所示。由圖6可知，A不變的條件下，Cd的振幅隨著f的增大而增大。當(dāng)脈動頻率f=8 Hz時，無量綱脈動振幅A分別為0.2、0.6、1時的Cd時間歷程曲線如圖7所示。由圖7可知，f不變的條件下，Cd的振幅隨著A的增大而增大。對比圖4、圖6可知，脈動流條件下圓柱阻力系數(shù)的振幅高出升力系數(shù)的振幅一個數(shù)量級，且阻力系數(shù)的變化速率也比升力系數(shù)的變化速率更快。與定常流條件下Cd的振幅[15]相比，脈動流條件下Cd的振幅更大，說明圓柱在脈動流條件下會受到更大的流向阻力，使圓柱體流向振動加強(qiáng)。

2.4 脈動頻率與脈動振幅對旋渦脫落頻率的影響

利用快速傅里葉變換(FFT)方法，可以將升力系數(shù)時間歷程曲線圖轉(zhuǎn)化為頻譜圖，使時域值轉(zhuǎn)化為頻域值，從而得到旋渦脫落頻率fv。A=0.6時，三種不同脈動頻率下(f=4、8、10 Hz)的升力系數(shù)Cl頻譜特性如圖8所示。圖8(a)中，當(dāng)f=4 Hz時旋渦脫落存在兩個主頻，主頻頻率分別為2.1、6.1 Hz；圖8(b)中，當(dāng)f=8 Hz時旋渦脫落存在三個主頻，主頻頻率分別為2.1、6.1、10.1 Hz；圖8(c)中，當(dāng)f=10 Hz時同樣存在三個主頻， 主頻頻率分別為2.1、8.1、12.1 Hz。由此可見，脈動流條件下圓柱繞流升力系數(shù)頻譜圖(圖8)中同時存在多個不同的主頻，其中振幅最大的主頻即為旋渦脫落頻率fv，且fv的振幅隨著脈動頻率f的增大而減小，其他主頻均為脈動頻率f與旋渦脫落頻率fv經(jīng)過相位疊加得到的頻率，這可能是由于渦的混合造成的。例如，在圖8(b)中，渦脫頻率為2.1 Hz，脈動頻率為8 Hz，兩個頻率經(jīng)過正負(fù)相位疊加分別得到了頻率為6.1、10.1 Hz的另外兩個主頻。

圖4 A=1時不同脈動頻率下的升力系數(shù)時間歷程曲線Fig.4 Time-history curve of lift coefficient at different pulsating frequency when A=1

圖5 f=10 Hz時不同無量綱脈動振幅下的升力系數(shù)時間歷程曲線Fig.5 Time-history curve of lift coefficient at different dimensionless pulsating amplitude when f=10 Hz

圖6 A=1時不同脈動頻率下的阻力系數(shù)時間歷程曲線Fig.6 Time-history curve of drag coefficient at different pulsating frequency (A=1)

圖7 f=8 Hz時不同無量綱脈動振幅下的阻力系數(shù)時間歷程曲線Fig.7 Time-history curve of drag coefficient at different dimensionless pulsating amplitude (f=8 Hz)

圖9 f=10 Hz時升力系數(shù)頻譜圖Fig.9 Frequency spectrogram of lift coefficient (f=10 Hz)

脈動頻率為f=10 Hz時，三種不同無量綱脈動振幅(A=0.2、0.6、1)的升力系數(shù)頻譜特性如圖9所示。圖9中存在3個主頻，其中振幅最大的即旋渦脫落頻率fv=2.1 Hz，另外兩個頻率是由脈動頻率f=10 Hz和旋渦脫落頻率fv=2.1 Hz相位疊加得到的，分別為8.1、12.1 Hz。由圖9還可以看出，頻譜圖中fv的振幅隨著無量綱脈動振幅A的增加而增加。

3 結(jié)論

通過研究脈動流條件下，脈動頻率和無量綱脈動振幅對圓柱繞流尾渦脫落、升阻力系數(shù)、旋渦脫落頻率等的影響規(guī)律，得出以下主要結(jié)論。

(1)當(dāng)脈動頻率增大時，尾渦長度變短，旋渦脫落速度加快；當(dāng)無量綱脈動振幅增大時，主導(dǎo)渦脫落速度減緩，另一側(cè)渦的脫落速度加快。

(2)隨著脈動頻率與無量綱脈動振幅的增大，升力系數(shù)與阻力系數(shù)的振幅均增大，且后者的振幅比前者的振幅高一個數(shù)量級，變化速率后者也更快。

(3)升力系數(shù)頻譜圖存在多個主頻，其中振幅最大的頻率為旋渦脫落頻率，其余頻率均是由脈動頻率與旋渦脫落頻率相位疊加的結(jié)果；隨著脈動頻率的增加，或無量綱脈動振幅的減小，旋渦脫落頻率的振幅減小。