周 強(qiáng)，曹曙陽，王 通，周志勇

（同濟(jì)大學(xué) 土木工程防災(zāi)國家重點(diǎn)實(shí)驗(yàn)室，上海 200092）

自然界的風(fēng)場主要有良態(tài)風(fēng)場和特異風(fēng)場。現(xiàn)行的橋梁與結(jié)構(gòu)抗風(fēng)設(shè)計(jì)方法主要是針對常規(guī)良態(tài)風(fēng)場，而實(shí)際風(fēng)災(zāi)害絕大多數(shù)是由極端風(fēng)氣候、突變風(fēng)等特異風(fēng)引起的，如龍卷風(fēng)、臺風(fēng)、下?lián)衾妆╋L(fēng)等。特異風(fēng)通常具有強(qiáng)速度剪切特性、特殊風(fēng)速剖面、強(qiáng)旋轉(zhuǎn)氣流以及復(fù)雜的時(shí)空分布特性。目前國內(nèi)針對特異風(fēng)的研究很少。瞿偉廉等［19］通過數(shù)值方法模擬了作用在某輸電塔結(jié)構(gòu)上的下?lián)舯┝骱奢d時(shí)程。趙楊等［18］利用主動(dòng)控制風(fēng)洞模擬出了雷暴沖擊風(fēng)陣風(fēng)剖面（如圖1所示）；然后，采用“階躍流法”模擬了風(fēng)速突變的時(shí)程；最后將高層結(jié)構(gòu)模型置于該突變氣流中，觀測在特殊氣流中結(jié)構(gòu)表面風(fēng)壓以及結(jié)構(gòu)空氣動(dòng)力學(xué)參數(shù)的變化特征。本文針對其中風(fēng)速剪切特性，對置于剪切流中的圓柱體氣動(dòng)力特性、尾流特性展開研究。

圖1 雷暴沖擊風(fēng)風(fēng)速剖面理論和試驗(yàn)值［18］Fig.1 The wind profile of thunderstorm and Boundary layer

目前針對考慮速度剪切效應(yīng)的圓柱體繞流問題的研究主要集中在航空航天領(lǐng)域，但由于其所關(guān)心的問題是升阻比，和風(fēng)工程存在差異，因此需要從風(fēng)工程角度對置于剪切流中的圓柱體繞流問題進(jìn)行研究。在這方面國外學(xué)者已有一些試驗(yàn)研究，Adachi和Kato［1］研究了在 Re=2.67×103～1.07×104和0＜β＜0.05范圍內(nèi)的圓柱體繞流問題，結(jié)果表明阻力和升力都隨著β的增加而增加，且升力的方向是由高速側(cè)指向低速側(cè)。Hayashi和 Yoshino［4］對在 Re=6 ×104，β=0.25 的流場中的圓柱氣動(dòng)力問題進(jìn)行了研究，發(fā)現(xiàn)阻力相比于均勻流是減少的，但升力方向仍然是由高速側(cè)指向低速側(cè)。Sumner和 Akosile［12］在小的剪切參數(shù)范圍（0.02 ＜β＜0.07）和雷諾數(shù) Re=4.0 ×104～9.0 ×104條件下的研究發(fā)現(xiàn)阻力與升力的結(jié)果和Hayashi和Yoshino［4］的結(jié)果基本相同。但是目前對這一問題的數(shù)值模擬研究較少，且都是進(jìn)行二維的計(jì)算。Tamura等［14］對處于Re=40，80，0 ＜β＜0.20 的剪切流中的圓柱體繞流問題進(jìn)行了數(shù)值模擬研究。Lei等［9］采用上游差分法對Re=50～160，0＜β＜0.2的流動(dòng)進(jìn)行了二維數(shù)值模擬。Lei和Tamura都預(yù)測在剪切流中前端的駐點(diǎn)位置會向高速側(cè)移動(dòng)；但是Lei的結(jié)果表明升力方向是指向低速側(cè)的，而Tamura則認(rèn)為是相反的方向。由此可見，對這種流場還有許多問題沒有澄清，因此很有必要采用三維數(shù)值模擬方法對其進(jìn)行進(jìn)一步的細(xì)節(jié)研究。

本文針對上述4個(gè)氣流特點(diǎn)中的風(fēng)速剪切特性，對置于剪切流中的圓柱體氣動(dòng)力特性展開研究，并且引入剪切參數(shù)β表示速度剪切強(qiáng)度，β=G（D/Uc），其中G是速度梯度，D是圓柱直徑，Uc是圓柱正前方來流的平均風(fēng)速，如圖2所示。雷諾數(shù)為Re=UcD/ν，因此當(dāng)Uc不變時(shí)，可以通過改變速度梯度G來改變剪切參數(shù)，即可以在雷諾數(shù)不變的情況下改變剪切參數(shù)。本文計(jì)算的雷諾數(shù)分別為 Re=60，80，150，200，220，500 和 1 000，其中當(dāng) Re=60，80，150，200和220時(shí)采用三維直接數(shù)值模擬（DNS），Re=500和1 000時(shí)采用基于Smagorinsky動(dòng)態(tài)亞格子模型的三維大渦模擬（LES）。本文主要研究氣動(dòng)力特性隨剪切參數(shù)β和雷諾數(shù)的變化情況，并分析產(chǎn)生這些氣動(dòng)現(xiàn)象的物理機(jī)理。

圖2 剪切流示意圖Fig.2 Schematic of shear flow configuration

1 控制方程與建模

1.1 控制方程

本文采用廣義曲線坐標(biāo)系來表示流體控制方程。DNS的連續(xù)性方程和N-S方程如式（1）和（2）所示：

式中：

式中：J為物理空間和計(jì)算空間之間Jacobian行列式，Um為廣義曲線坐標(biāo)下的逆變速度分量，Gmn為網(wǎng)格劃分的偏度張量。

LES的連續(xù)性和N-S方程為：

式中：～表示測試過濾操作。如Germano等［3］建議的，這里唯一的參數(shù)——測試過濾與格子過濾比率取2.0。同時(shí)為避免出現(xiàn)數(shù)值不穩(wěn)定，在計(jì)算中負(fù)的SGS渦粘性系數(shù)取為0。

1.2 網(wǎng)格劃分和算法

本文用Zang等［13］開發(fā)的基于 Co-located網(wǎng)格的數(shù)值模擬方法，在曲線坐標(biāo)系下，求解三維非定常不可壓縮N-S方程。在計(jì)算空間中，以控制體的中心作為定義點(diǎn)來定義速度分量和壓力，而以相應(yīng)的單元壁面中點(diǎn)為定義點(diǎn)來定義體積通量。用四階中心差分和四階中心插值離散對流項(xiàng)，用二階中心差分和二階中心插值離散擴(kuò)散項(xiàng)。關(guān)于時(shí)間發(fā)展，對流項(xiàng)采用Adams-Bashforth差分格式，而擴(kuò)散項(xiàng)則采用半隱式的Crank-Nicolson格式。

首先，笛卡爾坐標(biāo)系下中間速度分量可由式（11）得到：

通過坐標(biāo)轉(zhuǎn)換可以得到單元中心的逆變速度分量，從而插值得到單元壁面上的逆變速度分量。采用SOR法求解Poisson方程得到下一步的勢場：

然后，本文可以根據(jù)式（13）、式（14）得到笛卡爾坐標(biāo)系的速度分量和逆向速度分量：

上述算法已廣泛應(yīng)用于湍流場的數(shù)值模擬中。

1.3 數(shù)值模型和邊界條件

本文采用計(jì)算區(qū)域及邊界條件如圖3所示。

圖3 計(jì)算區(qū)域和邊界條件Fig.3 Computational domain and boundary conditions

圓柱體表面邊界條件：速度ui無滑移條件和擬壓φ的Neumann條件。

來流邊界條件：如圖2所示，采用速度剪切流作為來流，即U=1+Gy，v=0和w=0；擬壓的Neumann條件。同時(shí)為了避免在來流中出現(xiàn)反向流，當(dāng)y＜-1/G時(shí)，來流速度取為0。

展向邊界條件：速度和擬壓應(yīng)滿足周期性邊界條件。

2 數(shù)值驗(yàn)證

比較圓柱體近端和遠(yuǎn)端尾流的一階和二階湍流特性是理想的驗(yàn)證方法。但由于缺少在這雷諾數(shù)范圍內(nèi)的試驗(yàn)數(shù)據(jù)，因此本文首先在Re=50～1 000區(qū)域內(nèi)比較了Strouhal數(shù)和阻力系數(shù)，結(jié)果如圖4所示。然后圖5比較了本文在Re=1 000時(shí)得到的湍流特性結(jié)果、Kravchenko等［8］用大渦模擬得到的結(jié)果（Re=3 900）以及Ong和Wallace［11］等的試驗(yàn)結(jié)果。由此可見，本文進(jìn)行的直接模擬和大渦模擬是有效和可靠的。

3 數(shù)值結(jié)果及討論

本節(jié)主要研究了速度剪切流對圓柱繞流尾流結(jié)構(gòu)特征及氣動(dòng)力特性的影響。

3.1 繞流結(jié)構(gòu)特征

圖6給出了在不同的雷諾數(shù)下，Strouhal數(shù)隨剪切參數(shù)的變化情況?？梢钥闯?，Strouhal數(shù)基本不隨剪切參數(shù)變化，其大小只取決于雷諾數(shù)的大小。這一結(jié)果與 Sumner和 Akosile［12］、Cao 等［2］的試驗(yàn)結(jié)果相一致；與 Lei［9］、Kang［5］的數(shù)值模擬結(jié)果有細(xì)微差別，其原因可能在于他們進(jìn)行的是二維數(shù)值模擬。

圖7給出了瞬時(shí)尾流結(jié)構(gòu)隨剪切參數(shù)的變化情況，可以發(fā)現(xiàn)，高、低速兩側(cè)脫落的旋渦的不對性隨著剪切參數(shù)的增加而增大，Karman渦脫被抑制。

駐點(diǎn)漂移是剪切流中重要的現(xiàn)象之一，圖8給出了駐點(diǎn)位置隨著剪切參數(shù)變化的情況，其中θ0表示駐點(diǎn)位置所對應(yīng)的圓心角（高速側(cè)為正）?？梢钥闯觯S著β的增大駐點(diǎn)位置不斷向高速側(cè)移動(dòng)，且變化基本和β成線性關(guān)系。同時(shí)還可以從圖9看出，在低雷諾數(shù)下，駐點(diǎn)位置隨著雷諾數(shù)的增加而顯著改變，但是當(dāng)Re＞500后，其位置基本保持不變。所以本文認(rèn)為在較大雷諾數(shù)下，駐點(diǎn)位置基本不再受雷諾數(shù)效應(yīng)的影響，而主要由剪切系數(shù)β決定。

本文發(fā)現(xiàn)一個(gè)重要現(xiàn)象——兩側(cè)分離點(diǎn)隨剪切參數(shù)的變化規(guī)律。圖10給出了圓柱體兩側(cè)分離點(diǎn)位置隨剪切參數(shù)的變化情況?？梢钥闯龈咚賯?cè)的分離點(diǎn)向下游移動(dòng)，而低速側(cè)的分離點(diǎn)則向上游移動(dòng)，且隨著剪切參數(shù)呈線性變化。由于駐點(diǎn)和分離點(diǎn)位置的漂移，因此圓柱體周圍流場特征也發(fā)生了改變，分離邊界層的強(qiáng)度和厚度以及分離邊界層中產(chǎn)生的旋渦在圓柱兩側(cè)也不相同，從而在兩側(cè)形成不對稱的流動(dòng)。這些流場特征都將會影響圓柱體上氣動(dòng)力的產(chǎn)生，這些本文將在下面進(jìn)一步討論。

圖10 分離點(diǎn)位置隨剪切參數(shù)的變化圖（符號意義同圖6）Fig.10 Variation of separation point with shear parameter

3.2 氣動(dòng)力特性

作為一個(gè)例子，圖11給出了Re=1000時(shí)的圓柱體高速側(cè)和低速側(cè)的壓力系數(shù)分布圖?？梢钥闯觯瑑蓚?cè)的壓力分布不對稱，而且明顯可以發(fā)現(xiàn)駐點(diǎn)的位置向高速側(cè)移動(dòng)。同時(shí)還可以發(fā)現(xiàn)，壓力分布會隨著剪切參數(shù)變化。根據(jù)兩側(cè)壓力的相對大小，可見在從駐點(diǎn)到最小壓力點(diǎn)之間區(qū)域（A區(qū)），其特點(diǎn)為高速側(cè)的壓力更大；在由最小壓力點(diǎn)到分離點(diǎn)之間區(qū)域（B區(qū)），低速側(cè)的壓力更大；在兩側(cè)流動(dòng)的相互作用區(qū)域（C區(qū)），所以兩側(cè)壓力幾乎相同。A區(qū)的壓力分布主要是由駐點(diǎn)位置移動(dòng)引起的；而B區(qū)的壓力分布主要是由于整個(gè)流場中的速度剪切和分離點(diǎn)位置移動(dòng)的共同作用而產(chǎn)生的。在其他雷諾數(shù)下的壓力分布隨剪切系數(shù)的變化情況和Re=1000下的情況基本一致。

圖11 不同剪切參數(shù)下的圓柱體周圍壓力分布及其分區(qū)Fig.11 Variation of pressure distribution with shear parameter and its zoning

在剪切流中由于壓力分布不對稱性，導(dǎo)致升力產(chǎn)生。在B區(qū)，由于速度剪切的影響，低速側(cè)的壓力較大，從而產(chǎn)生一個(gè)由低速側(cè)向高速側(cè)作用的升力；而在A區(qū)，由于駐點(diǎn)漂移的影響，高速側(cè)的壓力較大，從而產(chǎn)生一個(gè)由高速側(cè)向低速側(cè)作用的升力。由于A區(qū)的作用更顯著，所以綜合的升力方向是由高速側(cè)指向低速側(cè)，如圖12 所示。Sumner等［12］、Cao 等［2］的試驗(yàn)研究，以及Lei等［9］數(shù)值模擬也都發(fā)現(xiàn)了類似的結(jié)論，但由于長寬比和紊流水平的不同，本文與上述文獻(xiàn)中的升力變化率有所不同。

圖12 升力系數(shù)隨剪切參數(shù)的變化Fig.12 Variation of lift coefficient with shear parameter at different Reynolds numbers

4 結(jié)論

本文采用三維直接數(shù)值模擬（DNS）和大渦模擬（LES）研究了雷諾數(shù)Re=60～1 000范圍內(nèi)速度剪切特性對圓柱繞流問題的影響，得到以下有益結(jié)論：

（1）在剪切流中的駐點(diǎn)向高速側(cè)移動(dòng)，且其對應(yīng)的圓心角隨剪切參數(shù)的增大而增大。

（2）分離點(diǎn)在高速側(cè)向下游移動(dòng)，在低速側(cè)向上游移動(dòng)，且分離點(diǎn)移動(dòng)的角度隨剪切參數(shù)的增加而增大，同時(shí)發(fā)現(xiàn)雷諾數(shù)對分離點(diǎn)移動(dòng)的影響很大。

（3）駐點(diǎn)以及分離點(diǎn)位置的移動(dòng)導(dǎo)致圓柱表面壓力發(fā)生變化，且這種變化對升力的貢獻(xiàn)要大于兩側(cè)速度不對稱性的貢獻(xiàn)，所以升力的作用方向由高速側(cè)指向低速側(cè)，而且隨著剪切參數(shù)的增加，升力系數(shù)也隨之增加。

（4）在Re=60～1 000范圍內(nèi)，Strouhal數(shù)基本不隨剪切參數(shù)變化，其大小取決于Re數(shù)的大小，而且隨著剪切參數(shù)的增加，Karman渦脫將會被抑制。

［1］Adachi T，Kato E.Study on the flow about a circular cylinder in shear flow［J］.Transactions of the Japan Society for Aeronautical and Space，1975，256：45-53.

［2］ Cao S，Ozono S，Hirano K，et al.Vortex shedding and aerodynamic forces on a circular cylinder in linear shear flow at subcritical Reynolds number［J］.Journal of Fluids and Structures，2007，23：703-714.

［3］Germano M，Piomelli U，Moin P，et al.A dynamic subgridscale viscosity model［J］.Physics of Fluids，1991，3（7）：1760-1765.

［4］Hayashi T，Yoshino F.On the evaluation of the aerodynamic forces acting on a circular cylinder in a uniform shear flow［J］.Transactions of the Japan Society for Mechanical Engineering，1990，56（552）：31-36.

［5］ Kang S.Uniform-shear flow over a circular cylinder at low Reynolds numbers［J］.Journal of Fluids and Structures，2006，22：541-555.

［6］Kim J，Moin P.Application of a fractional-step method to incompressible Navier-Stokesequations［J］. Journalof Computational Physics，1985，59：308-323.

［7］Kwon T S，Sung H J，Hyun J M.Experimental investigation of uniform shear flow past a circular cylinder［J］.ASME Journal of Fluids Engineering，1992，114：457-460.

［8］Kravchenko A G，Moin P.Numerical studies of flow over a circular cylinder at Re=3 900［J］.Physics of Fluids，2000，12（2）：403-417.

［9］Lei C，Cheng L，Kavanagh K.A finite difference solution of the shearflow overa circularcylinder［J］. Ocean Engineering，2000，27：271-290.

［10］Lilly D K.A proposed modification of the Germano subgridscale closure model［J］.Physics of Fluids，1992，4（4）：633-635.

［11］ Ong L，Wallace J.The velocity field of the turbulent very near wake of a circular cylinder［J］.Exp.Fluids，1996，20：441-453.

［12］Sumner D，Akosile O O.On uniform planar shear flow around a circular cylinder at subcritical Reynolds number［J］.Journal of Fluids and Structures，2003，18：441-454.

［13］ Zang Y，Street R L，Koseff J R.A non-staggered grid，fractional step method for time-dependent incompressible Navier-Stokes equationsin curvilinearcoordinates［J］.Journal of Computational Physics，1994，114：18-33.

［14］Tamura H，Kiya M，Arie M.Numerical study of viscous shear flow past a circular cylinder［J］.Transactions of the Japan Society for Mechanical Engineering，1994，46（404）：555-564.

［15］ Tamura T，Kitagishi T.Application of the interpolation method in generalized coordinate system to wake flows around a circular cylinder［J］.J.Struct.Constr.Eng.，AIJ，2001，545：27-34.

［16］ Williamson C H K.Vortex dynamics in the cylinder wake［J］.Annu.Rev.Fluid Mech，1996，28：477-539.

［17］ Zhang H Q，F(xiàn)ey U，Noack B R.On the transition of the cylinder wake［J］.Phys.Fluids，1995，7（4）：779-794.

［18］趙 楊，曹曙陽，Tamura Y，等.雷暴沖擊風(fēng)模擬及其的風(fēng)洞試驗(yàn)研究［J］.振動(dòng)與沖擊，2009，28（4）：1-3.

［19］瞿偉廉，王錦文.下?lián)舯┝鞯娘L(fēng)荷載模擬［J］.武漢理工大學(xué)學(xué)報(bào)，2008，30（2）：70-74.