雙向流固耦合的圓柱式橋墩繞流的三維數(shù)值模擬

任 輝，周亦唐，李 睿

(昆明理工大學(xué)建筑工程學(xué)院，云南昆明650000)

橋梁是道路跨越河流的主要形式，它起到保障公路運(yùn)輸暢通和宣泄洪水的作用［1］。橋墩建成后，橋位水域阻力增加，過水?dāng)嗝鏈p小，從而對周圍的水動(dòng)力環(huán)境產(chǎn)生一定影響。當(dāng)流體繞過物體時(shí)，在一定的情況下物體后部會(huì)發(fā)生漩渦脫落現(xiàn)象，特別是由此產(chǎn)生的周期性交替的作用力，會(huì)誘發(fā)結(jié)構(gòu)體振動(dòng)，導(dǎo)致疲勞破壞。為了提高橋梁抗水毀的能力，研究了橋墩與橋墩周圍三維繞流流場的相互耦合作用，對橋梁的安全具有非常重要的意義。

張兆順，等［2］指出結(jié)構(gòu)體繞流問題呈現(xiàn)出明顯的三維效應(yīng)。結(jié)構(gòu)體繞流不但涉及流動(dòng)分離、渦旋的生成和脫落，還包括旋渦的相互干擾等許多基本的問題，掌握流體流過結(jié)構(gòu)體障礙物的特性對于許多工程設(shè)計(jì)問題非常重要，而且對許多實(shí)際工程問題也有非常重要的指導(dǎo)作用。筆者以云南省漾濞(躍進(jìn))至云龍(諾鄧)二級(jí)公路1合同躍進(jìn)段K0+570大橋?yàn)檠芯康墓こ瘫尘埃ㄟ^系統(tǒng)的計(jì)算機(jī)仿真模擬，分析在一定水深、流速、流向等條件下橋墩周圍流場的分布情況，為建立橋梁抗水毀能力的評(píng)價(jià)體系和方法提供參考依據(jù)。

1 漩渦脫落的形成與機(jī)理

橋墩周圍的流場結(jié)構(gòu)主要包括墩前水面涌波、橋墩迎水面向下水流和尺度很大的漩渦。因圓柱表面存在邊界層，在柱體上發(fā)生的一個(gè)重要現(xiàn)象時(shí)邊界層分離(圖1)，邊界層外的流體可近似認(rèn)為是理想流體，因此沒有能量損失。邊界層外流體由于沒有動(dòng)能損失，仍在向前流動(dòng)，這樣在此區(qū)域形成流體的旋轉(zhuǎn)運(yùn)動(dòng)，亦即成漩渦，使流體不再貼著柱體表面流動(dòng)，而是從柱體表面分離出去，并形成向下游延展的自由剪切層。漩渦是柱體左右兩側(cè)交替地、周期性地發(fā)生的。當(dāng)一個(gè)漩渦向下游泄放，即自柱體脫落并向下游移動(dòng)時(shí)，它對柱體的影響及相應(yīng)的升力也隨之減小直到消失，而下一個(gè)漩渦又從對面一側(cè)發(fā)生，并產(chǎn)生同前一個(gè)相反方向的升力。因此，每一“對”漩渦具有相反的升力，并共同構(gòu)成一個(gè)垂直于流向的交變力的周期，引起結(jié)構(gòu)振動(dòng)［3-4］。

圖1 圓柱體表面邊界層分離Fig.1 Boundary layer separation of cylinder surface

2 數(shù)學(xué)模型

對不可壓縮黏性流體，在直角坐標(biāo)系下，其運(yùn)動(dòng)規(guī)律可用連續(xù)性方程和動(dòng)量方程來描述［5］。

1)連續(xù)性方程(質(zhì)量守恒方程)。在流場中，流體不斷的流進(jìn)或流出控制體，控制體中所包含的質(zhì)量可能隨時(shí)間變化。但是無論流動(dòng)如何復(fù)雜，就控制體而言，必須滿足質(zhì)量守恒定律，即單位時(shí)間內(nèi)通過控制面流入控制體的質(zhì)量之和等于單位時(shí)間內(nèi)控制體中質(zhì)量的增量：

式中：Vx，Vy，Vz分別為流體在3個(gè)方向上的速度，m/s。

2)運(yùn)動(dòng)方程(動(dòng)量守恒方程)。運(yùn)動(dòng)著的流體在任意時(shí)刻和任意空間點(diǎn)處的流體質(zhì)點(diǎn)在各外力作用下，將會(huì)以一定方式改變其運(yùn)動(dòng)學(xué)參數(shù)。將動(dòng)量守恒定律用于運(yùn)動(dòng)著的黏性流體質(zhì)點(diǎn)上，即可得到諸流動(dòng)參數(shù)間的特定關(guān)系，其數(shù)學(xué)表達(dá)形式即為Navier-Stokes方程(簡稱N-S方程)：

式中：ρ為流體密度，kg/m3;μ為流體動(dòng)力黏性系數(shù)，Pa·s;f為單位質(zhì)量流體體積力，N。

3 工程算例

3.1 模型建立及邊界條件

選用漾濞(躍進(jìn))—云龍(諾鄧)二級(jí)公路1合同躍進(jìn)段上一座40～20 m預(yù)應(yīng)力T形連續(xù)梁橋?yàn)槔?。?0孔采用結(jié)構(gòu)簡支橋面連續(xù)，后30孔采用結(jié)構(gòu)先簡支后結(jié)構(gòu)連續(xù)，全橋均為5孔一聯(lián)，共8聯(lián)，設(shè)9道伸縮縫;下部結(jié)構(gòu)，橋墩采用獨(dú)柱橋墩，橋臺(tái)采用重力式橋臺(tái)(圖2)。對其圓形獨(dú)柱式橋墩在三維的空間里進(jìn)行流體繞流與雙向流固耦合的數(shù)值模擬。

圖2 全橋立面展示(單位：cm)Fig.2 Elevation of the whole bridge

為更好地求解耦合方程，必須合理地定義好流場的初始條件和邊界條件。筆者參考已公開發(fā)表試驗(yàn)結(jié)果的圓柱模型［6-8］，定義如下。

1)結(jié)構(gòu)初始條件：在剛性柱壁面邊界條件為流固耦合面處X，Z方向無滑移，Y方向設(shè)為Free，并設(shè)置彈性約束條件E=300 N/m3(允許Y方向產(chǎn)生位移)。將圓柱表面與流體接觸的面設(shè)為流固耦合面，即流體計(jì)算結(jié)果與結(jié)構(gòu)計(jì)算結(jié)果相互耦合的聯(lián)系。

2)流體的初始條件：流場參考?jí)簭?qiáng)1 atm，水溫25℃。流場均勻來流速度為u=0.26 m/s。

3)遠(yuǎn)場條件：流場遠(yuǎn)處(流場出口)，平均靜態(tài)壓強(qiáng)P=0 Pa。流場兩壁面采用無滑移壁面。圓柱周圍的流場利用ANSYS CFX計(jì)算流體力學(xué)軟件求得，計(jì)算采用有限體積法及SIMPLE算法，離散格式采用Second Order Upwind二階迎風(fēng)格式。

3.2 模型尺寸及網(wǎng)格劃分

上游入口須取得足夠遠(yuǎn)以保證流場充分發(fā)展，下游出口離橋墩的距離也應(yīng)取得足夠大以保證流場達(dá)到穩(wěn)定［9］。參考已公開發(fā)表的模型試驗(yàn)確定模型尺寸如下所述，其中外部流域?yàn)?00 cm×100 cm，水深為20 cm，流速為0.26 m/s。橋墩直徑 D為10 cm。建模時(shí)流體區(qū)域用長方體模擬，橋墩用圓柱體模擬。橋墩布置于距上游5D處，出口邊界選在距離圓柱中心15D處。流體左右兩邊界距橋墩幾何中心(視橋墩為軸對稱結(jié)構(gòu))相等。物理模型和計(jì)算網(wǎng)格的劃分如圖3、圖4，網(wǎng)格劃分時(shí)對橋墩周圍和流體交界面附近網(wǎng)格進(jìn)行加密(采用ICEM進(jìn)行網(wǎng)格劃分，加密區(qū)域采用分塊實(shí)現(xiàn))，以保證計(jì)算精度。

圖3 計(jì)算的物理模型(單位：cm)Fig.3 Physical model obtained by the calculation

圖4 計(jì)算網(wǎng)格劃分Fig.4 Meshes division obtained by calculation

4 結(jié)果討論

模型數(shù)值計(jì)算是在CFX求解器下進(jìn)行，數(shù)值計(jì)算中時(shí)間步長取為0.05 s。圖5為由于圓柱橋墩受到流體作用，使得橋墩質(zhì)點(diǎn)位于z=H/2處質(zhì)點(diǎn)坐標(biāo)為(45，0，10)和 z=4H/5 處質(zhì)點(diǎn)坐標(biāo)為(45，0，16)的斷面位移時(shí)程曲線在22 s左右位移呈現(xiàn)周期性變化，圖中Ys表示圓柱在Y方向的位移。如圖6，在t=28及30 s時(shí)，流場運(yùn)動(dòng)是相同的，而t=29 s時(shí)的流場與t=28 s時(shí)的具有明顯的對稱性。即物體后部兩側(cè)產(chǎn)生了周期性漩渦脫落。

圖5 位移時(shí)程曲線Fig.5 The displacement time history curve

圖6 質(zhì)點(diǎn)位于z=4H/5斷面流場速度云圖Fig.6 Velocity nephograph when the particle in the plane of z=4H/5

圖7為Re=260時(shí)圓柱升力系數(shù)CL和阻力系數(shù)CD隨時(shí)間的變化曲線圖。文中，Z方向平均的升力和阻力系數(shù)分別定義為［10］：

式中：FL，F(xiàn)D分別為作用于圓柱上總的升力和阻力。

Strouhal數(shù)是描述圓柱繞流的一個(gè)重要參數(shù)，其定義為［10］：式中：St為Strouhal數(shù)，表征旋渦的脫離情況;fs為圓柱渦脫頻率;D為發(fā)生渦街柱體的迎風(fēng)面寬度(柱體迎面寬度即等于圓柱體直徑);U為來流速度。

圖7 Re=260時(shí)圓柱升力系數(shù)CL和阻力系數(shù)CD曲線Fig.7 Curve of lift coefficient CLand drag coefficient CDof circular cylinder at Re=260

由圖7阻力系數(shù)與升力系數(shù)時(shí)程曲線可以看出，升力系數(shù)時(shí)程曲線在32 s左右后顯現(xiàn)較為穩(wěn)定的周期性。對升力時(shí)程曲線進(jìn)行頻譜變換得到渦脫的主頻率為0.439 1 Hz，并由此計(jì)算得到Strouhal數(shù)為 0.168 8。

由圖8可見，在流動(dòng)折轉(zhuǎn)較大的區(qū)域，流動(dòng)速度及壓力增加，其最大值出現(xiàn)在橋墩迎水面兩側(cè)圓心角約70°～－90°范圍內(nèi)。由圖9可見在流體正對方向(即前駐點(diǎn))Cp為1，隨著向圓柱兩側(cè)擴(kuò)展，Cp值減小，在60°左右處達(dá)到最小，之后又快速增大到一個(gè)較穩(wěn)定值并在圓柱背后形成平坦的分布，由于流體力學(xué)將Cp開始平坦分布式的位置作為漩渦的脫離點(diǎn)，這與之前所述的邊界層分離即漩渦脫落現(xiàn)象較為吻合，其中Cp為平均壓力系數(shù);θ為圓柱表面測點(diǎn)所對應(yīng)的角度。

圖8 橋墩附近動(dòng)壓力分布Fig.8 Dynamic pressure distribution in the vicinity of pier

圖9 平均壓力系數(shù)周向分布Fig.9 Distribution of average pressure coefficients

5 結(jié)論

1)采用ANSYS CFX對橋墩附近流場進(jìn)行模擬，對Re=260的低雷諾數(shù)圓柱繞流，采用層流模型并結(jié)合合理的網(wǎng)格劃分，模擬圓柱表面周向壓力系數(shù)分布，預(yù)測分離點(diǎn)出現(xiàn)的位置，計(jì)算了阻力系數(shù)、升力系數(shù)及Strouhal數(shù)。數(shù)值計(jì)算的收斂曲線大概在700步(30 s)左右后顯示穩(wěn)定的周期性。

2)Re=260的低雷諾數(shù)圓柱繞流，由于圓柱表面渦脫在分離交界區(qū)的相互作用圓柱升力系數(shù)和阻力系數(shù)出現(xiàn)周期性的振動(dòng)。

［1］孫中堂.山區(qū)橋梁淺基礎(chǔ)水毀分析與防治［J］.山東交通科技，2002(3)：20-21.Sun Zhongtang.Analysis and prevention of shallow foundation of bridge in mountainous area［J］.Journal of Shandong Transportation Science and Technology，2002(3)：20-21.

［2］張兆順，崔桂香，許春曉.湍流理論與模型［M］.北京：清華大學(xué)出版社，2005.Zhang Zhaoshun，Cui Guixiang，Xu Chunxiao.Theory and Modeling of Turbulence［M］.Beijing：Tsinghua University Press，2005.

［3］張阿漫，戴紹仕.流固耦合動(dòng)力學(xué)［M］.北京：國防工業(yè)出版社，2002.Zhang A’man，Dai Shaoshi.Fluid Structure Interaction Dynamics［M］.Beijing：National Defense Industry Press，2010.

［4］高東光.公路與橋梁水毀防治技術(shù)［M］.北京：人民交通出版社，2002.Gao Dongguang.Innovative Technology for Highway and Bridge Flood Damage Prevention［M］.Beijing：China Communications Press，2002.

［5］江春波，張永良，丁則平.計(jì)算流體力學(xué)［M］.北京：中國電力出版社，2007.Jiang Chunbo，Zhang Yongliang，Ding Zeping.Computational Fluid Dynamics［M］.Beijing：China Electric Power Press，2007.

［6］Dong S，Karniadakis G E.DNS of flow past a stationary and oscillating cylinder at Re=10 000［J］.Journal of Fluids and Structures，2005，20(4)：519-531.

［7］Lam K，Wang F K，So R M C.Three-dimensional nature of vortices in the near wake of a wavy cylinder［J］.Journal of Fluids and Structures，2004，19(6)：815-833.

［8］Briaud J L，Chen H C，Wang J.Pier and Contraction Scour in Cohesive Soils［R］//National Cooperative Highway Project.Washington D C：Transportation Research Board，2004.

［9］凌建明，林小平，趙鴻鐸.圓柱形橋墩附近三維流場及河床局部沖刷分析［J］同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2007，35(5)：582-586.Ling Jianming，Lin Xiaoping，Zhao Hongduo.Analysis of three-dimensional flow field and local scour of riverbed around cylindrical pier［J］.Journal of Tongji University：Natural Science，2007，35(5)：582-586.

［10］白樺，李加武，夏勇.低雷諾數(shù)圓柱繞流數(shù)值模擬及控制措施［J］.建筑科學(xué)與工程學(xué)報(bào)，2010，27(4)：39-42.Bai Hua，Li Jiawu，Xia Yong.Numerical simulation and control measures of flow around circular cylinders at low Reynolds number［J］.Journal of Architecture and Civil Engineering，2010，27(4)：39-42.