不同倒角半徑四柱體繞流數(shù)值模擬及水動力特性分析

于定勇，康 驍，趙建豪

(中國海洋大學 工程學院，山東 青島 266100)

繞流是指流體繞過置于流體中的物體或是物體在流體中移動的過程，是生活中常見的一種物理現(xiàn)象。在實際工程中，柱體多以多柱體結(jié)構(gòu)形式布置，對倒角柱體而言，不同的柱體布置形式，柱體后方的流場變化、柱體受力情況也有很大不同。

不同倒角半徑柱體繞流的湍流場結(jié)構(gòu)復雜，國內(nèi)外已有學者開展了部分的研究工作。Tamura等[1]對不同倒角半徑柱體進行了數(shù)值模擬，發(fā)現(xiàn)柱體的角部形狀略微改變會極大影響柱體流場變化。Hu等[2]、Kumur等[3]使用粒子圖像流速測量技術(shù)(PIV)，研究了不同倒角半徑柱體的近尾流場特性和渦脫落周期的變化情況，發(fā)現(xiàn)前倒角半徑變化同后倒角相比對近尾流的影響要大。Tamura和Miyagi[4]對不同倒角半徑柱體進行了風洞模型試驗，發(fā)現(xiàn)相對于方柱，柱體截面形式為倒角時，平均阻力和尾流寬度有明顯減小。Zhang和Samtaney[5]通過直接數(shù)值模擬對不同倒角半徑柱體繞流進行了三維數(shù)值模擬，發(fā)現(xiàn)倒角半徑R+=0.125時，柱體的分離剪切層達到最長，升力系數(shù)均方根Cl-rms相較于方柱突降90%以上。

目前，四柱體布置下的柱體繞流研究工作主要針對圓形柱體及矩形柱體，對帶圓形倒角的柱體的研究較少。倒角柱體能減小應力集中[6]，降低結(jié)構(gòu)失效的可能性。Lam和Lo[7]在雷諾數(shù)Re=2 100方形四柱體布置條件下，分析了間距比和來流角度對柱體近尾流態(tài)和渦脫落頻率的影響。Lam和Fang[8]分析了不同間距比變化對方形四柱體的受力影響，根據(jù)結(jié)果總結(jié)了柱體阻力系數(shù)的變化規(guī)律。Lam等[9]研究了方形四圓柱在不同間距比下的渦脫落規(guī)律，發(fā)現(xiàn)圓柱剪切層的發(fā)展對渦脫落頻率的變化起著關(guān)鍵作用。殷長山等[10-11]通過數(shù)值模擬，對四圓柱在不同間距比、不同來流角度和不同雷諾數(shù)條件下的繞流流場特性進行了研究，研究發(fā)現(xiàn)間距比的變化會影響上游柱體剪切層的發(fā)展狀態(tài)，上游柱體剪切層的不穩(wěn)定性也會導致尾流形態(tài)的變化。

已有的研究工作表明，四柱體繞流的影響因素有很多，例如柱體間的排列方式、間距比、來流角度和雷諾數(shù)、倒角半徑等。以往對四柱體繞流的研究主要針對方柱和圓柱柱體進行，而對四柱體布置下的不同倒角半徑柱體繞流的研究較少。使用Fluent軟件，采用大渦模擬方法研究了在Re=3 900條件下，方形四柱體布置下(L=3D)6種不同倒角半徑三維柱體的后方瞬時流場、受力特性、時均流場的變化規(guī)律。

1 計算模型及參數(shù)設(shè)置

1.1 控制方程

假設(shè)流體為黏性不可壓縮，即流體密度不變，又假設(shè)溫度變化不大，因此能量方程可以忽略，故N-S方程只考慮連續(xù)方程和動量方程，引入濾波函數(shù)到N-S方程后可得：

(1)

(2)

1.2 模型參數(shù)

參考文獻[14]，選取計算模型和網(wǎng)格劃分如圖1所示，柱1、2為上游柱體，柱3、4為下游柱體；柱體直徑D=0.01 m，柱體間距L=3D，來流角度均為正向來流，順水流方向(x方向) 長度為40D，垂直水流方向(y方向) 長度為30D，計算域展向高(z方向)長為πD。不同倒角半徑柱體以R+=R/D來表示，其中R為倒圓角半徑，D為柱體邊長，設(shè)6種工況，分別為R+=0.0～0.5這6種截面形式，倒角半徑柱體截面形式如圖2所示。圖3表示四柱體布置下監(jiān)測線分布情況，在展向H=πD/2處設(shè)置截面，以前后柱體的中心處為x軸的起始位置，設(shè)置橫流方向監(jiān)測線6條，其中x/D=-0.5，0.0和0.5位于柱體之間，其他監(jiān)測線位于下游柱體后方，得到不同x/D處的順流向及橫流向速度分布。

圖1 計算區(qū)域圖及邊界層網(wǎng)格

圖2 不同倒角半徑柱體截面

圖3 方形四柱體布置下柱體監(jiān)測線布置

入口設(shè)置為速度入口，速度U0=0.391 m/s，水的運動黏度υ=1.003×10-6m2/s，得到雷諾數(shù)Re=U0D/υ=3 900。在每個柱體1.5D范圍內(nèi)使用Otype網(wǎng)格進行加密，考慮計算效率將展向節(jié)點數(shù)設(shè)為40個，首層厚度為0.003D，對應的Y+約為0.88，網(wǎng)格總數(shù)為3 734 192，計算區(qū)域總阻塞率為13.33%。設(shè)置時間步長取為 ΔtU0/D=0.039 1，以保證在計算域的庫朗數(shù)小于1。利用有限體積法離散控制方程，殘差精度設(shè)置為10-4。

1.3 邊界條件及初始條件

1) 進口條件：均勻來流速度(inlet)(U=U0，V=0，W=0)；

2) 出口條件：自由出流(outflow)；

3) 柱體表面：無滑移固壁面邊界(wall)(U=0，V=0，W=0)；

4) 上下左右邊界：對稱邊界(symmetry)。

2 模型驗證

為驗證所選取計算模型的可信性，下面對計算模型中的單柱體繞流和四柱體布置下的柱體繞流進行模型驗證。

圖4 單柱體繞流計算區(qū)域示意

2.1 單柱體模型驗證

單圓柱模型驗證選取的工況為Re=3 900下單圓柱繞流，圖4為單柱體繞流驗證計算區(qū)域圖。首先對近壁面網(wǎng)格和尾流區(qū)網(wǎng)格無關(guān)性進行分析，選用兩套網(wǎng)格進行考察，網(wǎng)格詳細設(shè)置如表1所示。表2給出了本文數(shù)值結(jié)果和已有單圓柱繞流試驗和物模結(jié)果的對比。由表2可知，這里兩套網(wǎng)格計算出的參數(shù)結(jié)果非常接近，得到的平均阻力系數(shù)Cd-ave、St數(shù)和回流長度同已有的PIV及數(shù)模結(jié)果吻合度較高，說明選取的網(wǎng)格和數(shù)值模型是可信的，考慮計算效率，文中所有網(wǎng)格設(shè)置與A1工況網(wǎng)格設(shè)置相同。

表1 網(wǎng)格設(shè)置情況

表2 數(shù)值計算結(jié)果驗證

2.2 四柱體模型驗證

圖5 平均阻力系數(shù)積分曲線

在驗證單柱體模型可信性之后，繼續(xù)對方形四柱體模型進行驗證。首先要確定四柱體模型合適的計算時間總長以及合理的采樣周期。為確定合適的計算時間總長，采用Zhang和Samtaney[5]中的平均阻力系數(shù)積分公式(式(3))，對R+=0.2～0.5時柱1的阻力系數(shù)進行積分，通過平均阻力系數(shù)積分曲線平穩(wěn)程度來衡量柱體后方繞流是否已經(jīng)穩(wěn)定。圖5為柱1的平均阻力系數(shù)積分曲線，可以看到隨著計算時長的增加，當tU0/D>400時曲線已經(jīng)基本水平穩(wěn)定，因此文中在tU0/D=450開始采樣進行時均統(tǒng)計。學者們選取的采樣周期數(shù)各不相同，Parnaudeau等[18]選取了420個渦脫落周期，Kravchenko等[17]采取了7個渦脫落周期，本文每隔0.001 s 統(tǒng)計一次，共采樣20～30個渦脫落周期。

(3)

四柱體模型驗證的結(jié)構(gòu)為方形布置四圓柱，雷諾數(shù)Re=1 500，間距L=3D，總阻塞率均為13.33%，驗證對比工況為殷長山[11]的PIV物理模型試驗。圖6表示在x/D=1.0、1.5和2.0處順流向時均速度和殷長山[11]PIV試驗下游不同橫斷面處速度分布情況對比。

圖6 下游不同橫斷面處時均順流速度分布(Re=1 500)

由圖6可得在三個不同橫斷面處，計算的速度曲線同PIV試驗測點速度走勢吻合較好，說明本文的四柱體模型對柱體近尾流場模擬較為準確，上文設(shè)置的模型和流場是可信的。因此在此模型基礎(chǔ)上，分別對6種不同倒角半徑柱體在方形四柱體布設(shè)條件下進行了三維數(shù)值模擬，并分析了其時均和瞬時流場特性、漩渦脫落形態(tài)和受力特性的變化規(guī)律。

3 模擬結(jié)果及分析

3.1 柱體瞬時流場特性分析

本節(jié)研究了方形四柱體布置情況下，柱體后方的瞬時流場變化及近尾流結(jié)構(gòu)隨倒角半徑改變的變化規(guī)律。 圖7表示tU0/D=475時四柱體布置下的瞬時渦量圖，可知在R+=0.0時，四柱體之間的渦流變化較為混亂，在上游柱體轉(zhuǎn)角處剪切層進行初始分離，并在上游柱體后方破碎并部分附著于下游柱體表面，后方柱體間的尾流相互交雜并一同向下游繼續(xù)發(fā)展。

圖7 四柱體展向截面瞬時渦量等值線圖(1/s)

圖8 四柱體三維渦量側(cè)視圖(Q=10)

在R+=0.1時，發(fā)現(xiàn)剪切層在上游柱體分離后在柱體間形成渦街，下游柱體后方尾流寬度明顯增加(圖8所示)；由于橫向間柱體干擾較強，后方渦街成對稱分布，這種渦街對稱分布在R+=0.5時也有體現(xiàn)。在R+=0.2～0.4時，從上游柱體分離的剪切層得到發(fā)展，并附著于下游柱體前端及兩側(cè)，此時剪切層應為“再附著模式”，尾流寬度較于R+=0.1時明顯減小，且可以看到后方渦街沒有成對稱分布，皆呈現(xiàn)同側(cè)分布。為了更直觀表現(xiàn)出柱體瞬時三維渦量，采用Hunt等[20]建議的Q準則來表示，圖中Q值均為10。從圖8可以看到隨著R+的增加，上游柱體剪切層初始分離處速度值開始逐漸減小，剪切層的破碎分離均在上游柱體后下游柱體前，其中在R+=0.0、0.1時剪切層破碎較早，使得上下游柱體間產(chǎn)生渦街導致尾流寬度的增加。

3.2 柱體受力特性分析

本節(jié)研究了當柱體倒角半徑變化時，四柱體布置情況下柱體升阻力系數(shù)的變化規(guī)律。因柱1和柱2、柱3和柱4的受力情況基本相同，選取柱1和柱4來反映柱體的受力情況。如圖9所示，研究發(fā)現(xiàn)同單柱體繞流類似[6]，隨著R+的增大上游柱體的Cd-ave逐漸減小，且在R+=0.5時相較于R+=0.0時Cd-ave減小約50%以上；下游柱體Cd-ave總體要小于上游柱體，除了在R+=0.1處產(chǎn)生突增并達到正值以外，其余均隨R+增加而緩慢增大。由圖10可以看出相較于上游柱體，下游柱體的Cl-rms總體偏大；其中上、下游柱體的Cl-rms在R+=0.1時均產(chǎn)生突增，原因可能是R+=0.1時剪切層破碎過早，在上、下游柱體之間產(chǎn)生渦街，導致了柱體橫向受力較大。

圖9 柱體平均阻力系數(shù)分布

圖10 柱體均方根升力系數(shù)分布

圖11 柱體St數(shù)分布圖

在不同倒角半徑下，四個柱體的St數(shù)都相等，與殷長山[11]得出的L/D>1.5時四柱的渦脫落頻率在同一工況下基本一致的結(jié)論相吻合。圖11以柱1為例，給出柱體St數(shù)隨R+的變化規(guī)律。可以看到當R+=0.0四柱體均為方柱的時候St數(shù)最小，隨著R+的增大各柱體St數(shù)的逐漸增大，R+=0.5時St數(shù)小幅降低。當間距比都為L/D=3時，并列方柱兩柱的St數(shù)均為0.142，串列方柱兩柱的St數(shù)均為0.127，方形布置四方柱的St數(shù)均為0.141，可以看到隨著R+的增大，方形四柱體的渦脫落頻率比并列方柱和串列方柱渦的脫落頻率要更高。

圖12和圖13分別表示不同倒角半徑下各柱體的升力系數(shù)和阻力系數(shù)時程曲線。從圖12可以明顯發(fā)現(xiàn)，在R+=0.0和0.2～0.5時，上游柱體的升力系數(shù)曲線波動幅度顯著小于下游柱體；在R+=0.1、0.5時，上游兩柱體柱1柱2的升力系數(shù)曲線相位相反，同時下游兩柱體柱3和柱4的升力系數(shù)曲線相位也相反。觀察圖7可以發(fā)現(xiàn)，R+=0.1、0.5時柱體后方渦街成對稱分布；在R+=0.2、0.3和0.4時，上游兩柱體下游兩柱體的升力系數(shù)曲線相位均相同，從圖7可以發(fā)現(xiàn)此時柱體后方的渦脫落傾向獨立性發(fā)展，分布模式為同向分布；而在R+=0.0時，四個柱體的升力系數(shù)曲線相位之間沒有明顯的對應關(guān)系，既不是對稱分布也不是同向分布，可能原因是此時柱體近尾流附近渦脫落比較不規(guī)律。從圖13阻力系數(shù)時程曲線可以看出，上游柱體阻力系數(shù)明顯要大于下游柱體較且變化幅度很小，上游、下游阻力曲線皆沒有明顯的周期性;上游柱體的阻力系數(shù)隨著R+的增大逐漸減小，下游柱體的波動區(qū)間值除了在R+=0.0時突然增加，其余皆隨著R+的增大逐漸增大，這與圖9的內(nèi)容一致。

圖12 四柱體布置下柱體升力系數(shù)時程曲線

圖13 四柱體布置下柱體阻力系數(shù)時程曲線

3.3 柱體時均流場特性分析

本節(jié)主要研究四柱體布置情況下倒角半徑變化時，對柱體后方時均流場的影響。從圖14可以看到，R+=0.1柱體的順流向速度波動相較于其他柱體明顯不同，在其它R+下的柱體順流向波動區(qū)域主要集中于下游柱體前端及后方，且影響范圍均較小，而R+=0.1時則分布在各個柱體的后方且波動區(qū)域較大，其原因參考圖4可以發(fā)現(xiàn)在R+=0.1時各個柱體后方均有渦街形成，這使得各柱體后方速度變化較為劇烈。

從圖15可以看出，上游柱體的前方和下游柱體后方的回流區(qū)內(nèi)的順流向速度波動不明顯，表明此處順流向速度變化較為穩(wěn)定；另外可以發(fā)現(xiàn)，順流向速度的最小值除了R+=0.1時分布在上游柱體后方，其他R+下多分布于下游柱體前側(cè)的外端，這在圖16中也可觀測到。

從圖16(a)、(b)、(c)看出，在上下游柱體之間，順流向速度為“雙V型”分布，其中除了R+=0.1的柱體，其他柱體順流向速度最低點均向外側(cè)偏移，結(jié)合圖15可以發(fā)現(xiàn)，此時順流向速度最低處多集中于下游柱體前端的外側(cè)，而在R+=0.1時則主要分布在上游柱體正后方。在x/D≥3以后，各倒角半徑柱體速度分布趨向于一致，均從“深雙V型”向“淺雙V型”過渡。

圖14 順流向均方根速度分布云圖(z=πD/2)

圖15 順流向時均速度分布云圖 (z=πD/2)

圖16 不同x/D條件下順流向時均速度

圖17表示橫流向均方根速度分布情況，可以看到在除了在R+=0.1時橫流向速度波動程度都比較小，產(chǎn)生原因是四柱體間的相互作用限制了尾流在橫流向的發(fā)展；另一方面由于上游柱體倒角的影響，使得柱間剪切層得到增加，從而抑制了柱體間渦街的產(chǎn)生。另外需要注意的是，在R+=0.1時發(fā)現(xiàn)四個柱體后方橫流向速度波動較大，參照圖7(b)可以看到此時上下游柱體后方均發(fā)生渦街脫落，再一次說明渦街脫落是影響橫流向速度波動的主因。

從圖18橫流向時均速度分布云圖可以看到，由于來流作用于柱體前端并向兩側(cè)分流，使得上游柱體前端橫向流速分布一正一負，皆是上端為正下端為負；除了R+=0.1時，其他情況下下游柱體后端橫向流速分布也為一正一負，與上游柱體前端相反，皆是上端為負下端為正。

圖17 橫流向均方根速度分布云圖(z=πD/2)

圖18 橫流向時均速度分布云圖(z=πD/2)

圖19表示不同x/D條件下橫流向時均速度分布，在上下游柱體之間，即在x/D=-0.5，0.0和0.5處，橫流向速度波動多集中在y/D=±3之間且為反對稱分布。其中在x/D=-0.5和0.0時，R+=0.1柱體的橫流向速度分布相對其他R+有較大區(qū)別，結(jié)合圖18發(fā)現(xiàn)在上游柱體后方出現(xiàn)一對反向的速度波動區(qū)，這是由于在此處產(chǎn)生渦街所導致，而渦街的產(chǎn)生會使橫流向速度值較大并發(fā)生劇烈波動。從下游柱體后方的監(jiān)測線(x/D=3.0～10.0)可以看到，各倒角半徑柱體的橫流向速度分布規(guī)律大致相同，說明柱體截面形式的變化對后方流場影響較小，隨著x/D的增大，波動趨于平緩，并在x/D=10.0處趨于水平。

圖19 不同x/D條件下橫流向時均速度

4 結(jié) 語

利用Fluent軟件采用大渦模擬方法對Re=3 900時方形四柱體布置條件下(L=3D)6種不同倒角半徑柱體繞流過程進行了研究，分析了數(shù)值模擬結(jié)果，發(fā)現(xiàn)倒角半徑的變化對柱體后方流場特性及柱體受力有明顯的影響。

1)R+=0.1時剪切層在上游柱體分離后并在柱體間形成渦街，而其他R+柱體剪切層均再附著于下游柱體；在R+=0.1時，由于剪切層分離后在上、下游柱體間形成渦街，導致下游柱體后方尾流寬度明顯增加。

2) 隨著R+增大，上游柱體的Cd-ave逐漸減小，下游柱體的Cd-ave除了在R+=0.1處產(chǎn)生突增以外，均隨R+增大而增加；下游柱體的Cl-rms總體大于上游柱體，上、下游體的Cl-rms走勢相同，均在R+=0.1時產(chǎn)生突增。

3)R+=0.1、0.5時，上下游兩柱體的升力系數(shù)時程曲線相位相反，此時柱體后方的渦脫落為對稱分布，而在R+=0.2、0.3和0.4時，上下游兩柱體的升力系數(shù)時程曲線相位相同，此時下游柱體后方的渦脫落傾向獨立性發(fā)展。

4)R+=0.1柱體的順流向速度波動分布在各個柱體的后方且波動區(qū)域較大，在其它R+下的柱體順流向波動區(qū)域主要集中于下游柱體前端及后方，且影響范圍均較?。怀薘+=0.1時各柱體橫流向速度波動程度很小，對比R+=0.1柱體發(fā)現(xiàn)渦街脫落是影響橫流向速度波動的主因。