錯(cuò)列圓柱的三維大渦模擬

汪 洋,陳 立

(四川省公路規(guī)劃勘察設(shè)計(jì)研究院有限公司，四川成都 610000)

斜拉索、電纜等鈍體結(jié)構(gòu)由于其柔度大、阻尼小易導(dǎo)致各種風(fēng)致振動(dòng)問題，對(duì)于這類問題一般將結(jié)構(gòu)簡(jiǎn)化為圓柱模型進(jìn)行研究[1-5]。圓柱繞流是流體力學(xué)中的經(jīng)典問題之一，當(dāng)粘性流體流經(jīng)鈍體時(shí)，會(huì)發(fā)生邊界層分離現(xiàn)象，進(jìn)而在尾部形成周期性渦脫，同時(shí)伴隨著鈍體兩側(cè)壁面壓力的變化，從而對(duì)結(jié)構(gòu)產(chǎn)生有規(guī)律的氣動(dòng)力作用。特別對(duì)于雙圓柱繞流的情況，圓柱間的耦合干擾作用會(huì)導(dǎo)致其氣動(dòng)力特性有別于單圓柱繞流，因此雙圓柱繞流問題一直是專家學(xué)者研究的熱點(diǎn)[6-10]。

大部分研究集中在串列圓柱繞流上。Zdravkovich認(rèn)為串列圓柱之間的流體流動(dòng)形態(tài)與柱間中心距有關(guān)，當(dāng)兩圓柱中心距為1.0D～1.5D時(shí)，上游圓柱分離出的自由剪切層會(huì)越過下游圓柱；當(dāng)兩圓柱中心距為1.5D～4.0D時(shí)，剪切層會(huì)附著在下游柱體上；當(dāng)中心距大于4D時(shí)，上游和下游柱體會(huì)分別形成渦脫，其中D為圓柱直徑。Kitagawa運(yùn)用大渦模擬研究了Re=22000時(shí)不同中心距下的串列圓柱，結(jié)果表明：當(dāng)中心距小于3.25D時(shí)，只有下游圓柱出現(xiàn)漩渦脫落；當(dāng)中心距大于3.25D時(shí)，上下游均發(fā)生渦脫現(xiàn)象，并且發(fā)現(xiàn)中心距在3.25D附近時(shí)，上下游的平均阻力系數(shù)和升力系數(shù)會(huì)發(fā)生跳躍性變化。Tokoro與Komatsu的理論研究與風(fēng)洞試驗(yàn)表明:當(dāng)串列圓柱中心距為4.3D風(fēng)向角為15 °時(shí)，下游圓柱會(huì)出現(xiàn)明顯振動(dòng)。Brika和Lanevill采用風(fēng)洞試驗(yàn)的方法對(duì)中心距為10D～25D的串列雙圓柱進(jìn)行了遠(yuǎn)距失穩(wěn)現(xiàn)象的研究。也有少部分學(xué)者對(duì)并列圓柱的繞流特性開展了研究，Sumner[15]分別對(duì)小中心距(1.0D～1.2D)、中等中心距(1.2D～2.2D)、大中心距(>2.2D)的并列圓柱繞流進(jìn)行了分析。在實(shí)際工程中，鈍體的布置形式相對(duì)于流體流動(dòng)方向一般都是錯(cuò)列的，然而錯(cuò)列雙圓柱的繞流狀態(tài)卻很少受到研究者的關(guān)注。

本文采用大渦模擬的方法對(duì)Re=3900中心距3D的錯(cuò)列雙圓柱繞流進(jìn)行數(shù)值模擬，得到了不同風(fēng)向角下的升阻力系數(shù)、尾流相對(duì)速度以及風(fēng)壓系數(shù)等參數(shù)，并分析了下游圓柱升力系數(shù)出現(xiàn)極大值的原因。

1 錯(cuò)列雙圓柱繞流數(shù)值模型

1.1 計(jì)算工況

風(fēng)向角α的范圍為0～24 °，角度從0 °開始每2 °設(shè)置一個(gè)工況，一共13個(gè)工況，圓柱相對(duì)中心距W=3.0D(圖1)。

圖1 圓柱布置示意

1.2 計(jì)算域與邊界條件

計(jì)算域設(shè)置為：流向長(zhǎng)度20D，橫向長(zhǎng)度10D，展向長(zhǎng)度3.2D，以保證湍流尾流的充分發(fā)展。計(jì)算區(qū)域劃分為三個(gè)部分，均采用結(jié)構(gòu)化網(wǎng)格：第一部分為3D×3D×3.2D矩形近模型區(qū)域，該區(qū)域?yàn)榧用芫W(wǎng)格區(qū)，越靠近圓柱網(wǎng)格的密度越大，第一層網(wǎng)格到壁面的距離為2×10-4m，Y+值約為0.8；第二部分為直徑9D的圓形區(qū)域，方便在調(diào)整工況時(shí)改變風(fēng)向角；第三部分為外圍計(jì)算域，適當(dāng)降低網(wǎng)格密度以節(jié)省計(jì)算機(jī)資源。形成的網(wǎng)格見圖2，模型的總網(wǎng)格數(shù)量達(dá)到182×104個(gè)。

圖2 計(jì)算域網(wǎng)格劃分

計(jì)算域左側(cè)為速度入口，流速為U0=0.57m/s；右側(cè)為壓力出口；上下前后流域邊界為對(duì)稱邊界；圓柱表面為無(wú)滑移壁面。

坐標(biāo)系選擇：雙圓柱中心連線的中點(diǎn)作為坐標(biāo)原點(diǎn)，x軸方向與來(lái)流方向一致，y軸垂直于計(jì)算域上下邊界，z軸由右手定則確定。

1.3 參數(shù)設(shè)置

流體為空氣，密度ρ=1.225 kg/m3，動(dòng)力粘度μ=1.7894×10-5pa·s。亞格子尺度模型采用WALE模型，壓力速度耦合求解算法采用對(duì)大多數(shù)問題都適用的SIMLPE算法，梯度離散采用節(jié)省計(jì)算量的最小二乘法，壓力離散格式采用標(biāo)準(zhǔn)格式，動(dòng)量離散格式采用具有高階格式的有界中心差分格式。瞬態(tài)求解時(shí)時(shí)間步長(zhǎng)Δt=0.005s，以保證庫(kù)朗數(shù)小于5且每個(gè)時(shí)間步迭代收斂。

1.4 網(wǎng)格獨(dú)立性驗(yàn)證

為驗(yàn)證網(wǎng)格獨(dú)立性，對(duì)雷諾數(shù)Re=3900的單圓柱在同樣的網(wǎng)格參數(shù)(計(jì)算域尺寸、網(wǎng)格密度、時(shí)間步長(zhǎng)與錯(cuò)列圓柱模擬相同)條件下進(jìn)行了數(shù)值模擬，并與Lourenco[17]與Norberg[18]的風(fēng)洞試驗(yàn)進(jìn)行對(duì)比。從圖3與圖4中可見時(shí)均速度與表面壓力系數(shù)的曲線與試驗(yàn)曲線吻合較好。

圖3 單柱平均流速

圖4 單柱壓力系數(shù)

2 模擬結(jié)果

2.1 尾流平均風(fēng)速

在中心平面上(z=0)本文選取了3個(gè)剖面來(lái)監(jiān)測(cè)尾流平均風(fēng)速，分別是x/D=0，x/D=2.1，x/D=3。因工況較多，本文僅列出風(fēng)向角α=0°與α=24°情況下的模擬結(jié)果。

圖5給出了α=0°時(shí)平均流速沿y軸的分布規(guī)律。從圖中可以看出速度關(guān)于x軸對(duì)稱，這與結(jié)構(gòu)的對(duì)稱性相一致，隨著湍流的發(fā)展，速度分布曲線的形狀由“U”形逐漸變?yōu)椤癡”形，說明湍流尾流在周邊流體的影響下，速度逐漸趨近于周邊流體的速度。在y=0處，速度大小出現(xiàn)極小值，風(fēng)速大小從圓柱中心向兩側(cè)迅速增長(zhǎng)，快速增長(zhǎng)區(qū)域大約為2D，之后風(fēng)速趨于穩(wěn)定。圖6為α=0°平均流速云圖，可以直觀地看到下游圓柱完全浸沒在上游圓柱的尾流之中，氣流在圓柱之間形成了較大范圍的回流區(qū)，在下游圓柱尾部也形成了長(zhǎng)度約為1D的回流區(qū)。

圖5 α=0°平均流速分布

圖6 α=0°平均流速云圖

圖7給出了α=24°時(shí)平均流速沿y軸的分布規(guī)律。由于圓柱對(duì)氣流的阻擋作用，在y方向上與圓柱位置對(duì)應(yīng)處的流速會(huì)出現(xiàn)兩個(gè)極小值，下游圓柱尾流流速的極值出現(xiàn)在圓柱的正后方，而上游圓柱尾流因?yàn)殚g隙流的擠壓效應(yīng)，速度極值的出現(xiàn)位置會(huì)向y軸負(fù)向偏移。圖8為α=24°的平均流速云圖，下游圓柱不再浸沒在上游圓柱的尾流之中，圓柱之間出現(xiàn)間隙流并分別在尾部出現(xiàn)回流區(qū)，下游圓柱回流區(qū)長(zhǎng)度約為1.5D，上游圓柱回流區(qū)略小于上游圓柱，長(zhǎng)度約為1.3D。

圖7 α=24°平均流速云圖

圖8 α=24°平均流速云圖

2.2 平均風(fēng)壓系數(shù)

沿圓柱周向布置36個(gè)監(jiān)測(cè)點(diǎn)，監(jiān)測(cè)點(diǎn)間隔為10 °。圖9為α=0°時(shí)的圓柱時(shí)均壓力系數(shù)分布。從圖中可以發(fā)現(xiàn)上游圓柱表面正對(duì)來(lái)流處(θ1=0°)的壓力系數(shù)最大，大約為1.0。隨著氣流向圓柱兩側(cè)擴(kuò)展，Cp值迅速減小，并逐漸轉(zhuǎn)變?yōu)樨?fù)壓，在θ1=70°處達(dá)到最大負(fù)壓，大約為-1.0，說明氣流邊界層在該處與上游圓柱分離，分離角θsep=70°，邊界層分離后壓力系數(shù)Cp逐漸回升并增大到一個(gè)較為穩(wěn)定的值，在圓柱背后形成穩(wěn)定的負(fù)壓分布；下游圓柱受到柱間回流區(qū)域的影響，表面壓力系數(shù)均為負(fù)值，θ2=30°與θ2=330°處的負(fù)壓值最大，最大負(fù)壓值約為-0.77，在θ2=70°與θ2=290°處負(fù)壓值最小，最小值為-0.14，說明氣流從上游圓柱分離之后在該處進(jìn)行了再附，隨著氣流向圓柱后方延伸，負(fù)壓值略微增大之后趨于穩(wěn)定，并在下游圓柱后方形成穩(wěn)定的壓力分布區(qū)域。

圖9 α=0°平均風(fēng)壓系數(shù)

圖10為α=24°時(shí)的圓柱時(shí)均壓力分布。在正對(duì)來(lái)流處(θ=0°)兩圓柱的壓力系數(shù)最大，最大值均為0.25。隨著氣流在圓柱壁面上的流動(dòng)，壓力系數(shù)迅速減小，并逐漸轉(zhuǎn)變?yōu)樨?fù)壓，上游圓柱的氣流分離角為70 °與290 °，在該角度處的壓力系數(shù)均為-0.21，壓力系數(shù)在圓柱截面上下對(duì)稱分布；下游圓柱的氣流分離角為70 °與300 °，對(duì)應(yīng)的壓力系數(shù)分別為-0.20、-0.24，由于兩圓柱之間間隙流的影響，壓力系數(shù)在下游圓柱截面的分布是非對(duì)稱的。

圖10 α=24°平均風(fēng)壓系數(shù)

2.3 升、阻力系數(shù)與Strouhal數(shù)

圖11到圖14分別給出了α=0°與α=24°的升、阻力系數(shù)時(shí)程曲線。α=0°時(shí)力系數(shù)隨時(shí)間的變化不平穩(wěn)，但仍存在一定周期性特征。上游圓柱升力系數(shù)約為0，阻力系數(shù)約為0.93；下游圓柱升力系數(shù)約為0，阻力系數(shù)約為-0.24。由于上游圓柱尾流的撞擊作用，下游圓柱力系數(shù)的波動(dòng)幅值要大于上游圓柱。

圖11 α=0°上游圓柱力系數(shù)

圖12 α=0°下游圓柱力系數(shù)

圖13 α=24°上游圓柱力系數(shù)

圖14 α=24°下游圓柱力系數(shù)

α=24°時(shí)力系數(shù)隨時(shí)間的變化具有較強(qiáng)的周期性。上游圓柱升力系數(shù)約為-0.02，阻力系數(shù)約為1.05；下游圓柱升力系數(shù)約為-0.38，阻力系數(shù)約為0.98。

圖15展示了時(shí)均升力系數(shù)Cl與風(fēng)向角α之間的關(guān)系，可以發(fā)現(xiàn)風(fēng)向角的變化不會(huì)影響上游圓柱升力系數(shù)，其值穩(wěn)定在0附近。下游圓柱的升力系數(shù)均小于零，說明下游圓柱有向上游圓柱靠近的趨勢(shì)，升力系數(shù)曲線形狀近似于“V”字形，風(fēng)向角在8 °以下時(shí)，升力系數(shù)絕對(duì)值隨著風(fēng)向角的增加而增大，即下游圓柱向上游圓柱靠近的趨勢(shì)愈加明顯，當(dāng)風(fēng)向角大于8 °時(shí)，升力系數(shù)絕對(duì)值又隨著風(fēng)向角的增加而減小。

圖15 升力系數(shù)與風(fēng)向角的關(guān)系

圖16展示了時(shí)均阻力系數(shù)Cd與風(fēng)向角α之間的關(guān)系，上游圓柱阻力系數(shù)隨著風(fēng)向角的增大而緩慢增加。下游圓柱的阻力系數(shù)隨著風(fēng)向角的增加而明顯增大，阻力系數(shù)由負(fù)值轉(zhuǎn)變?yōu)檎担f明下游圓柱的運(yùn)動(dòng)趨勢(shì)由向上游運(yùn)動(dòng)逐漸轉(zhuǎn)變?yōu)橄蛳掠芜\(yùn)動(dòng)。

圖16 阻力系數(shù)與風(fēng)向角的關(guān)系

將升力系數(shù)時(shí)程進(jìn)行傅里葉變換并進(jìn)行無(wú)量綱化處理可以得到圓柱的Strouhal數(shù)。圖17給出了Strouhal數(shù)與角度α的關(guān)系，可以發(fā)現(xiàn)上下游圓柱的Strouhal數(shù)始終保持一致，α=0°時(shí)Strouhal數(shù)最小，最小值為0.169，對(duì)應(yīng)的渦脫頻率0.963 Hz；α=6°時(shí)Strouhal數(shù)最大，最大值為0.348，對(duì)應(yīng)的渦脫頻率1.985 Hz。

圖17 St數(shù)與風(fēng)向角的關(guān)系

2.4 渦脫形式

Q判據(jù)[19-20]為識(shí)別旋渦的一種方法，形式簡(jiǎn)單，并且運(yùn)用廣泛，所以本文以Q判據(jù)作為漩渦出現(xiàn)的判斷準(zhǔn)則，并給出了α=0°與α=24°時(shí)圓柱三維渦結(jié)構(gòu)。當(dāng)α=0°時(shí)，下游圓柱會(huì)阻礙上游尾渦的發(fā)展，導(dǎo)致兩圓柱之間形成許多小型漩渦，并且多個(gè)方向的漩渦交錯(cuò)在一起使得流態(tài)十分復(fù)雜見圖18；當(dāng)α=24°時(shí)，上游圓柱的尾渦會(huì)繞過下游圓柱，并與下游圓柱尾渦相結(jié)合，結(jié)合后的尾流渦結(jié)構(gòu)比α=0°的更寬且更加有規(guī)律性見圖19，這也解釋了24 °風(fēng)向角比0 °風(fēng)向角的的升阻力系數(shù)更具有規(guī)律性的原因。

圖18 α=0°三維渦結(jié)構(gòu)

圖19 α=24°三維渦結(jié)構(gòu)

2.5 大升力原因分析

觀察圖15可以發(fā)現(xiàn)，下游圓柱在α=8°時(shí)受到了很大的升力，為了進(jìn)一步探討下游圓柱受到很大升力的機(jī)理，選取α=8°的某個(gè)時(shí)刻來(lái)分析局部流場(chǎng)特性，見圖20中右側(cè)圓柱。圖21為α=8°時(shí)的流線圖/風(fēng)速云圖，圖22為α=8°時(shí)的風(fēng)壓云圖。

圖20 α=8°渦量云圖

圖21 α=8°流線圖/風(fēng)速云圖

圖22 α=8°風(fēng)壓云圖

由圖21可看出，靠近下游圓柱下表面的地方，出現(xiàn)了高速間隙流，最大速度系數(shù)達(dá)到了1.75，對(duì)應(yīng)圖22中相同的地方出現(xiàn)了最大負(fù)壓區(qū)域，該負(fù)壓區(qū)域使圓柱受到較大的升力作用，這是出現(xiàn)較大升力的第一個(gè)原因；其次，觀察圖21中圓柱下表面的流線，形成了一個(gè)分離泡，說明在間隙流的影響下，圓柱下側(cè)邊界層在此分離，分離后會(huì)形成多個(gè)小型漩渦，小型漩渦與分離泡會(huì)在圓柱表面造成局部強(qiáng)負(fù)壓，這就是導(dǎo)致圓柱出現(xiàn)較大升力的第二個(gè)原因；此外，觀察圖21圓柱上表面的流線，也存在邊界層分離現(xiàn)象，在圖22中對(duì)應(yīng)位置出現(xiàn)了負(fù)壓極值，分離角為79 °，根據(jù)之前的分析得知α=0 °時(shí)分離角為70 °，比較可知α=8°時(shí)負(fù)壓極值點(diǎn)上移，導(dǎo)致下游圓柱會(huì)受到更大的升力作用。

以上為下游圓柱在α=8°時(shí)受到很大的升力的三個(gè)原因，并總結(jié)如下：

(1)高速間隙流引起的強(qiáng)負(fù)壓。

(2)小型漩渦與分離泡產(chǎn)生的局部強(qiáng)負(fù)壓。

(3)邊界層分離點(diǎn)上移。

3 結(jié)論

在Re=3900的情況下采用大渦模擬的方法對(duì)中心距為3D、風(fēng)向角0～24 °的錯(cuò)列雙圓柱進(jìn)行了數(shù)值模擬，研究了空氣動(dòng)力特性隨風(fēng)向角的變化，分析了大升力出現(xiàn)的原因。主要結(jié)論如下：

(1)隨著風(fēng)向角的改變，錯(cuò)列圓柱會(huì)出現(xiàn)兩種繞流形態(tài)：0～8 °時(shí)，圓柱之間出現(xiàn)大范圍的回流區(qū)域，8～24 °時(shí)，柱間回流區(qū)域轉(zhuǎn)變?yōu)楦咚匍g隙流，這種轉(zhuǎn)變對(duì)圓柱氣動(dòng)力有顯著影響。

(2)隨著風(fēng)向角的增大，升、阻力曲線會(huì)變得更有規(guī)律性，且下游圓柱升、阻力曲線波動(dòng)幅值始終大于上游圓柱，這是因?yàn)樵诖箫L(fēng)向角時(shí)，上游圓柱尾渦會(huì)繞過下游圓柱并與下游圓柱尾渦相結(jié)合，形成更寬更規(guī)則的渦結(jié)構(gòu)。

(3)在8 °風(fēng)向時(shí)，下游圓柱升力系數(shù)會(huì)出現(xiàn)一個(gè)峰值，分析其原因主要有三點(diǎn)：一是高速間隙流引起的強(qiáng)負(fù)壓；二是小型漩渦與分離泡產(chǎn)生局部強(qiáng)負(fù)壓；三是邊界層分離點(diǎn)由70 °上升到79 °。