多圓柱繞流的數(shù)值模擬及力學(xué)分析

趙紅曉

（同濟(jì)大學(xué)航空航天與力學(xué)學(xué)院，上海 200092）

0 引言

圓柱的繞流現(xiàn)象在工程中廣泛存在，如管排［1-2］、橋墩［3-4］、電纜和高聳建筑［5］等結(jié)構(gòu)周圍的流動(dòng)問題。1940 年11 月7 日，美國華盛頓州塔克馬海峽大橋因風(fēng)振致毀，事故原因正是空氣繞流形成的卡門渦街引起的橋梁共振。作為典型鈍體的一種，近代學(xué)者對圓柱繞流的研究一直沒有停息過。

1911 年，馮·卡門觀察到在一定條件下的定常來流繞過某些物體時(shí)，物體兩側(cè)會(huì)周期性地脫落出旋轉(zhuǎn)方向相反、并排列成有規(guī)則的雙列線渦，并對渦街的誘導(dǎo)速度、穩(wěn)定性和阻力等做了分析［6］。Coutanceau等［7］揭示了圓柱后尾流形態(tài)和雷諾數(shù)之間的對應(yīng)關(guān)系。Zdravkovich［8］研究了雙圓柱串列在流體流動(dòng)速度很小時(shí)，不同的間距比L/D（圓柱間距與圓柱直徑之比）對兩柱間隙中的速度分布和下游柱尾流場中的速度分布具有明顯影響。Le 等［9］發(fā)現(xiàn)雙圓柱繞流在間距比L/D=2 時(shí)尾流產(chǎn)生的渦街迅速消失，間距比為4.5 時(shí)兩條尾流的形狀隨時(shí)間隨機(jī)發(fā)生變化，出現(xiàn)渦脫落方向相同、相反等4 種典型的尾流形狀。Yokoi等［10］研究了雙柱并列（非等直徑），間距比為0.25 ～1.5 之間時(shí)流體形態(tài)發(fā)生交替偏斜，間距比大于1.5時(shí)，兩柱后分別形成各自的渦脫落。張愛社等［11］利用有限元法對等邊三角形布置的等直徑3 圓柱繞流問題進(jìn)行了研究，表明小圓柱的位置、大小、個(gè)數(shù)以及雷諾數(shù)變化等對主圓柱有影響。Kuo 等［12］研究了2 個(gè)并列小圓柱置入尾流區(qū)，可有效地抑制渦街的產(chǎn)生，并使主圓柱的升力系數(shù)降低了70% ～80%，壓力、阻力減小也很明顯。張洪軍等［13］研究了附屬小圓柱對并列雙圓柱尾流流場的影響。Muddada 等［14］研究了放置于主圓柱后的旋轉(zhuǎn)小圓柱可以有效地抑制渦街的生成。

由于流動(dòng)的復(fù)雜性，目前以多圓柱群為對象的研究相對較少，工程情況中圓柱以圓柱群的形式出現(xiàn)的案例卻十分常見，本文選取不等直徑多圓柱群實(shí)驗(yàn)裝置進(jìn)行仿真研究。

1 數(shù)值模擬

1.1 控制方程

本文采用Fluent軟件進(jìn)行數(shù)值模擬，選用Simplec算法（半隱式方法Semi-Implicit Method for Pressure Linked Equations）進(jìn)行計(jì)算。本文研究中，低雷諾數(shù)下水的流動(dòng)，可近似為不可壓縮流體，則二維連續(xù)性方程可簡化為：

動(dòng)量方程為動(dòng)量守恒定律在流體力學(xué)中的表示形式：

式中：u、v、w為3 個(gè)速度分量；p為壓強(qiáng)；ρ為流體密度；μ為流體動(dòng)力黏度。

升力系數(shù)Cl和阻力系數(shù)Cd是流體力學(xué)中的無因次量，用來表示物體在流體中的升力和阻力大小

式中：Fl、Fd分別為物體受到的升力和阻力；v∞為無限遠(yuǎn)場的來流速度；S為物體在來流方向上的投影面積。

雷諾數(shù)

式中，v、L分別為流場的特征速度和特征長度。雷諾數(shù)物理上表示慣性力和黏性力量級的比。

本文均設(shè)置左側(cè)邊界為速度入口，右側(cè)邊界為壓力出口，遠(yuǎn)場壁面為對稱邊界條件，圓柱壁面采用無滑移壁面。流體介質(zhì)為水，密度ρ=1000 kg/m3，動(dòng)力黏度μ=1.0 ×10-3kg/（m·s）。

1.2 計(jì)算模型

多圓柱繞流研究實(shí)驗(yàn)裝置如圖1 所示，以還原實(shí)驗(yàn)裝置中的流動(dòng)規(guī)律，劃分模型進(jìn)行多組仿真計(jì)算。

圖1 流態(tài)演示裝置斷面示意圖

為避免過于龐大的模型造成不必要的計(jì)算時(shí)間浪費(fèi)，需對模型進(jìn)行適當(dāng)?shù)暮喕蛣澐帧?/p>

如圖2 所示，模型保留了實(shí)驗(yàn)裝置中近壁面繞流和多圓柱群繞流的流域，并在入口和出口處做了適當(dāng)?shù)募娱L。與實(shí)驗(yàn)室現(xiàn)象進(jìn)行對比，觀察能否正確模擬實(shí)際流動(dòng)，保留了實(shí)驗(yàn)裝置中的柵格，多圓柱排列方式、尺寸和編號如圖3 所示。

圖2 實(shí)驗(yàn)裝置物理模型及尺寸

圖3 多圓柱群尺寸及編號示意

1.3 網(wǎng)格劃分及計(jì)算方法驗(yàn)證

采用結(jié)構(gòu)化網(wǎng)格的方法進(jìn)行網(wǎng)格劃分，建立oblock與圓柱面以及柵格前緣、尾緣的圓角關(guān)聯(lián)，將其余邊與壁面和柵格長邊關(guān)聯(lián)，并在湍流情況中，運(yùn)用式（7）估算第1 層網(wǎng)格到柱面的距離Δy進(jìn)行相應(yīng)加密。最終得到的網(wǎng)格具有較好的質(zhì)量和長寬比，網(wǎng)格劃分如圖4 所示。

圖4 實(shí)驗(yàn)裝置網(wǎng)格劃分

網(wǎng)格劃分與層流驗(yàn)證模型時(shí)相同。為保證圓柱面上第1 層網(wǎng)格滿足y+≈1，估算確定第1 層網(wǎng)格到柱面的距離［15］

代入特征長度L=D=0.01 m，Re=3900，得到Δy=0.03 mm，加密后網(wǎng)格如圖5 所示。

圖5 圓柱周圍及局部網(wǎng)格加密圖

2 層流模型和紊流模型驗(yàn)證

參考文獻(xiàn)［16］，使用相同的物理模型進(jìn)行層流驗(yàn)證。左側(cè)為速度入口，右側(cè)為壓力出口，遠(yuǎn)場壁面為對稱邊界條件，圓柱壁面采用無滑移壁面。計(jì)算區(qū)域20D×25D，采用結(jié)構(gòu)化網(wǎng)格，計(jì)算參數(shù)見表1，計(jì)算結(jié)果見表2。

表1 層流驗(yàn)證模型仿真參數(shù)

表2 層流模型驗(yàn)證結(jié)果與文獻(xiàn)數(shù)據(jù)對比

由表2 可知，計(jì)算結(jié)果與文獻(xiàn)［15］基本一致，可使用層流模型對小雷諾數(shù)的情況進(jìn)行模擬。

參考文獻(xiàn)［15］中使用相同的物理模型進(jìn)行紊流驗(yàn)證。左側(cè)為速度入口，右側(cè)為壓力出口，遠(yuǎn)場壁面為對稱邊界條件，圓柱壁面采用無滑移壁面。計(jì)算區(qū)域20D×40D，靠近圓柱加密網(wǎng)格，采用結(jié)構(gòu)化網(wǎng)格，計(jì)算參數(shù)見表3，計(jì)算結(jié)果見表4。

表3 湍流驗(yàn)證模型仿真參數(shù)

表4 湍流模型驗(yàn)證結(jié)果與文獻(xiàn)數(shù)據(jù)對比

計(jì)算結(jié)果與文獻(xiàn)［15］基本一致，可使用分離渦模型進(jìn)行模擬。

3 仿真計(jì)算

分別采用以上驗(yàn)證計(jì)算模型和劃分網(wǎng)格的方法分別對多圓柱繞流模型（見圖3（a））進(jìn)行數(shù)值模擬，表5中雷諾數(shù)均為小圓柱雷諾數(shù)，參數(shù)設(shè)置見表5。

表5 實(shí)驗(yàn)裝置中多圓柱繞流仿真參數(shù)

由速度云圖圖6（a）可見，在Re=200 時(shí)，以大圓柱的前緣為中心，水流在駐點(diǎn)受阻后，被迫向兩邊流動(dòng)，在靠近圓柱尾部的邊界上，主流開始與圓柱體分離，邊界層分離后，在分離區(qū)的下游形成回流區(qū)，即為尾渦區(qū)。后流經(jīng)柵格區(qū)域，在流經(jīng)柵格的各狹縫中，由于截面積的減少，流動(dòng)皆有一定的加速。流體流經(jīng)多列小圓柱時(shí)，列與列之間的尾渦出現(xiàn)摻混而紊動(dòng)更強(qiáng)烈，渦的長度由于受到限制變短。第4 列圓柱則直接在后方產(chǎn)交替脫落渦。

由圖6（b）速度云圖可見，當(dāng)Re=5000 時(shí)，大圓柱尾渦出現(xiàn)明顯的兩列渦，且脫落形成長度明顯變短，渦街的寬度變窄，尾流呈復(fù)雜的湍流特征渦。小圓柱的第1 列圓柱限制了后面第2、3 列圓柱分離剪切層的卷起，使其不能產(chǎn)生強(qiáng)烈的旋渦脫落，第3、4 列圓柱的渦在后方混合，形成更為復(fù)雜的尾流。各圓柱后方仍有一定范圍的回流泡，旋渦脫落現(xiàn)象在下游第4 列圓柱后尾渦脫落成偏斜的尾流模式，隨著雷諾數(shù)的增大，第4 列后的紊動(dòng)就越強(qiáng)烈。

當(dāng)Re=10000 時(shí)，從圖6（c）的速度云圖可觀察到其流動(dòng)與Re=5000 時(shí)十分相似，但紊動(dòng)更強(qiáng)烈，尾流模式和圓柱受力不具有周期性，是由于尾渦區(qū)的長度和紊動(dòng)強(qiáng)度與來流的雷諾數(shù)有關(guān)，雷諾數(shù)越大，紊動(dòng)越強(qiáng)烈。

當(dāng)Re=200 時(shí)，由圖7（a）、（b）可見，大圓柱受到最大阻力系數(shù)最大，小圓柱中第1 列（從左向右看）圓柱阻力次之，第3 列圓柱所受阻力略小于第1 列圓柱，位于中心線上的圓柱12、32 仍相比同列受到更大的阻力。第2 和第4 列圓柱所受阻力比第3 列略小，第4列圓柱所受阻力略大于第2 列。第4 列圓柱升阻力系數(shù)的振幅比前幾列大。

圖7 多圓柱繞流阻力系數(shù)隨時(shí)間變化曲線

當(dāng)Re=5000 時(shí)，由圖7（c）、（d）和（e）可見，大圓柱的最大阻力系數(shù)仍是最大，但比Re=200 時(shí)小。圓柱12 阻力系數(shù)次之，第4 列圓柱的阻力系數(shù)比第3 列大，位于中心線上的圓柱12、32 仍相比同列其他圓柱受到更大的阻力，各圓柱的阻力系數(shù)均有一定幅度的振蕩，后方的圓柱較前方有更大的振蕩幅度。

當(dāng)Re=10000 時(shí)，由圖7（c）、（d）和（e）可見，各圓柱的阻力系數(shù)均有一定幅度的下降和振蕩，后方的圓柱較前方有更大的振蕩幅度。同樣位于中心線上的圓柱12、32 仍相比同列受到更大的阻力。

由圖7 中不同雷諾數(shù)下的阻力系數(shù)可見，大圓柱的升力系隨著雷諾數(shù)的增加而降低，小圓柱的繞流阻力系數(shù)總體上也隨著雷諾數(shù)的增加而降低，隨著雷諾數(shù)的增加，流動(dòng)由層流變?yōu)槲闪鳎膭?dòng)更強(qiáng)烈，各阻力系數(shù)均隨時(shí)間產(chǎn)生微小振蕩，振幅逐列增加。

當(dāng)Re=200 時(shí)，由圖8（a）～（d）可見，第4 列小圓柱升力系數(shù)大約在±4.5 這一范圍內(nèi)振蕩，其次是圓柱31、33，圓柱11、13。圓柱21、22 升力系數(shù)在零上下振蕩，對稱軸上的大圓柱和小兩圓柱12、32 所受升力為零。第1、3 列和第2、4 列同側(cè)圓柱所受升力仍具反向性。

圖8 多圓柱繞流升力系數(shù)隨時(shí)間變化曲線

當(dāng)Re=5000 時(shí)，由圖8（e）～（h）可見，升力系數(shù)的振蕩幅度也逐列增加，第4 列圓柱的升力系數(shù)振幅較大，最大值為其他圓柱的2 倍以上，第3 列升力系數(shù)次之，第1 列最小。不過前2 列圓柱的振蕩較小，除去中心線上的圓柱12、32，其他圓柱受力方向也沒發(fā)生改變，而大圓柱和小圓柱第4 列圓柱隨著尾流的復(fù)雜化，其升力系數(shù)開始較大幅度的振蕩，且方向會(huì)隨時(shí)間改變，說明紊動(dòng)強(qiáng)度增加，能量損失更大。

當(dāng)Re=10000 時(shí)，由圖8（i）～（l）可見，規(guī)律和雷諾數(shù)5000 時(shí)相似，升力系數(shù)的振蕩幅度較之雷諾數(shù)5000 時(shí)更大，流體紊動(dòng)強(qiáng)度繼續(xù)增加。

對比Re=5000 時(shí)的速度云圖和實(shí)驗(yàn)裝置的流動(dòng)圖像（見圖9），可發(fā)現(xiàn)兩者較好符合，認(rèn)為仿真結(jié)果與實(shí)際流動(dòng)是很接近的，可用數(shù)值模擬的方法替代實(shí)驗(yàn)裝置來演示紊流的流態(tài)。

圖9 實(shí)驗(yàn)裝置流動(dòng)圖像對比

4 結(jié)語

本文對實(shí)驗(yàn)裝置簡化模型在3 種雷諾數(shù)下進(jìn)行了數(shù)值計(jì)算和流態(tài)、阻力系數(shù)和升力系數(shù)分析，可知，多圓柱繞流在雷諾數(shù)Re=5000 時(shí)的速度云圖和實(shí)驗(yàn)裝置的流動(dòng)圖像符合，認(rèn)為仿真結(jié)果與實(shí)際流動(dòng)是很接近的，可采用模擬的方法去演示流體的紊動(dòng)機(jī)理，為流體力學(xué)虛擬仿真開發(fā)提供有益的探索。