基于有限體積法的群樁繞流數值模擬

莊 茜，李紹武，祁澤鵬

基于有限體積法的群樁繞流數值模擬

莊 茜1,2,3，李紹武1，祁澤鵬1

(1. 天津大學水利工程仿真與安全國家重點實驗室，天津 300072；2. 天津城建大學土木工程學院，天津 300384；3. 長沙理工大學水沙科學與水災害防治湖南省重點實驗室，長沙 410114)

采用COASTALTOOL軟件中的潮流模塊TC2D模擬樁繞流問題．通過細網格模型，研究樁斷面形狀、排列方式及樁心距對樁周圍水流的影響．對于粗網格模型中群樁阻力的影響，提出基于阻力疊加原理的有樁單元水流阻力計算公式；對于單元中群樁的遮流作用，提出遮流面積折合水深的概化方法．這種概化計算模式中考慮了群樁中各樁間的遮蔽影響．

樁繞流；群樁阻力；遮流；概化

由群樁組成的透空式建筑物，對周圍水流和泥沙運動的擾動要比實體墻建筑物小，在某些特定區(qū)域用做調整流場的工程措施往往具有一定優(yōu)勢，技術上需要根據群樁的阻流特性確定群樁布置方式和樁間距，理論上歸結為孤立建筑物繞流問題．

物體繞流問題十分復雜，至今只能對簡單、規(guī)則斷面形狀的圓樁或方樁在低雷諾數情況下給出解析解或半經驗理論解．對于群樁阻流特性的問題需借助物理模型試驗和數值計算．Zdravkovich[1]通過水槽試驗，研究了兩圓樁串列和交錯放置，且樁心距小于5倍柱徑時的繞流流場，結果表明，樁心距在3倍樁徑是臨界值，小于此值時，上游圓樁的渦尚未脫落，即進入第2個樁的影響范圍；而超過此臨界值時，上游樁后交替脫落的渦對下游樁的繞流力產生影響. Williamson[2]對兩并列圓樁的繞流尾跡進行了研究，通過流場顯示技術得到兩圓樁的尾跡變化．李景銀等[3]在低雷諾數下，在水槽中對4個正方形排列、樁間距為4倍樁徑的圓樁進行了繞流試驗，分別采用激光誘導熒光技術(LIF)和粒子成像測速(PIV)法測量了流場分布．鄧紹云等[4-5]通過水槽試驗，探討了規(guī)則排列群樁的柱阻力系數隨樁間距的變化．

在有網格數值模擬計算中，存在兩種基本方法．一種是按單樁解析尺度進行網格剖分，進行原樁數值計算；另一種是采用較大網格，對單元內多個樁進行概化后計算．原樁計算的可以獲得較高精度，但計算量大，概化計算的關鍵是如何對各樁的阻流效果進行合理概化．黃筱云等[6]采用自適應四叉樹網格下的N-S方程數值求解圓樁繞流，其結果與前人試驗結果一致．Stansby[7]用渦方法研究了雙圓樁并排、串列、交錯放置時的阻力特性，認為排列方式對阻力有影響．吳文權[8]通過數值計算研究了對稱來流條件下多圓樁繞流問題，認為雙圓樁繞流流動是雙穩(wěn)態(tài)的，多圓樁的排列方式對流場影響很大，前排對后排的影響尤為顯著．Farrant等[9]用數值模擬方法研究了雷諾數Re＝200時4個圓樁排列成正方形時的繞流問題，兩樁間距分別為3倍和5倍樁徑時，圓樁后方均出現旋渦脫落現象．張愛社等[10]用數值方法模擬了3圓樁布置成等邊三角形時的二維繞流問題，得出各圓樁阻力系數和流場分布，結果表明，間距較小時3圓樁相互干擾，流動有偏向下游某個圓樁的傾向．唐士芳等[11-12]提出用折合底部摩阻力的經驗方法來計算群樁繞流阻力，其中對于樁繞流阻力的處理缺乏明確的物理機制．解鳴曉等[13]提出在保證阻力等效的前提下以大尺度樁代替群樁的方法，概化中未考慮樁間影響．

本文首先通過細網格模型，研究不同樁斷面、群樁排列方式及樁心距對樁周圍水流的影響．對于粗網格模型中群樁的影響，提出一種綜合考慮底部摩阻力與樁繞流阻力的理論公式和樁間遮流效果的概化方法，并通過對比細網格原樁模擬和粗網格概化模擬結果，驗證該概化方法的有效性．

1 二維水流數學模型

運用自行開發(fā)的海岸軟件COASTALTOOL中的二維潮流模塊TC2D進行水流數值計算，模型的基本方程及求解方法概述如下．

1.1 控制方程

在如圖1所示的坐標系下，淺水波運動方程為式中：H為全水深，H＝η－Zb，η為水位，Zb為底高程；u、v分別為x、y方向上的流速分量；g為重力加速度；τx、τy分別為水流阻力在x、y方向的分量；ρ為水的密度；f為Coriolis頻率參數，為地球自轉頻率，為當地緯度；νh為水平方向的紊動黏滯系數，采用Smagorinsky亞格子紊動模型計算，即

式中：Cs為系數，取0.1～0.2；A為單元面積．

圖1 坐標系示意Fig.1 Sketch of coordinate system

1.2 邊界條件

開邊界處給水位過程．水平方向固邊界采用“不穿透”條件，沿固邊界的法向分量恒為零，即

式中：v為流速矢量；n為固邊界的單位外法線矢量．

1.3 控制方程的數值求解方法

采用三角形和四邊形混合單元，以提高計算效率．離散變量定義在各單元的中心，單元交界處變量產生間斷，交界面數值通量采用Osher格式計算，有關方程離散過程詳見文獻[14]．

1.4 群樁概化方法

將群樁所在單元視為水單元，采用等效阻力系數和等效過水面積相結合的方法對群樁進行概化，具體方法如下．

1.4.1 單元綜合水流阻力

單元內有樁時，單元綜合水流阻力可視為底部摩阻力與樁繞流阻力之和，即

式中：Ai為第i個單元的面積；wiT和piT分別為底部摩阻力和樁繞流阻力．底部摩阻力wiT的表達式為

式中：Awi為扣除樁位后單元i的面積；C為謝才系數，可用曼寧公式計算；上標n表示時間層．

樁繞流阻力piT為單元中各樁繞流阻力之和，即

式中：Mi為單元i內樁的個數；αfm為考慮群樁內各樁相互影響后的繞流阻力校正系數，根據文獻[5]和筆者的大量數值研究結果，αfm與樁的排列方式、樁心距等因素有關，其值取決于各樁對綜合阻力的貢獻，可根據試驗結果或數值試驗結果選取；Cd為單樁繞流阻力系數；Ap為單樁在水流方向投影的面積，為單樁在水流方向投影的寬度為第i單元的自然水深．

假設單元內樁尺寸統(tǒng)一，離散方程中的x向單元綜合水流阻力分量可表示為

y方向單元綜合水流阻力分量可表示為

1.4.2 過水面積等效

群樁不但增加了水流阻力，也改變了局部過水面積．過水面積等效是將單元中樁群遮擋的過水面積折合成垂向高度，在單元水深中予以扣除．需扣除的折合水深計算公式為式中：αdm為考慮樁間遮蔽作用的調整系數；Lm為第i單元中第m根樁在單元中沿水流方向影響域的長度(見圖2)；Bm為第i個單元中第m根樁在單元中沿水流方向的影響域的寬度．

遮蔽調整系數αdm的取值目前尚無理論方法．從工程實用角度，初步可按如下方法取值．若單元內沿水流方向單樁或并列多根樁時αdm取1；若橫向多根樁不并列(錯位布置)，則對位于上游的首個樁αdm取1，其后的樁則根據文獻[1, 9-10]中有關前樁對后樁的影響僅限于其后3倍樁徑D范圍的結論，按式αdm=min(δ/3D,1)計算，其中δ 為前、后樁在水流方向投影距離．往復流可按不同方向算出的折合水深中最不利方向考慮．

圖2 樁在單元中沿水流方向影響域示意Fig.2Sketch of influence domain of pile in an element along flow direction

2 模型有效性驗證

采用Roulund等[15]的水槽試驗資料對二維水流數學模型的有效性進行檢驗．試驗為單個圓樁繞流問題．數值水槽長1,000,m，寬200,m，縱比降為1.17/10,000，底部糙率與物模一致，取0.022．水深為0.54,m，圓樁直徑為0.536,m，置于水槽橫斷面中心，均勻流流速為0.326,m/s，樁雷諾數為1.75×105．采用三角形網格，樁附近進行網格加密(見圖3)．計算區(qū)域左端為來流邊界，右端為出流邊界，上、下側均為滑移固邊界．

樁前、后垂向平均流速模擬結果與實測結果對比如圖4所示，樁前模擬結果與實測結果符合很好，樁后實測值呈現先增大、再減小、而后又增大的過程，這可能是由于樁后1倍D范圍內為尾流區(qū)，渦動使得流場變化復雜，故二維模型無法準確模擬樁后較近范圍內的復雜流場，但在尾流區(qū)以外的計算結果與實測結果是吻合的．

圖3 圓樁繞流計算區(qū)域及網格劃分Fig.3 Computational domain and mesh discretization of flow around cylindrical pile

根據伯努利方程導出樁前水位壅高/Z DΔ的計算公式為

式中Fr為弗勞德數．

仍采用上述數值水槽，計算條件如表1所示．駐點A處(見圖3)相對水位壅高/Z DΔ隨弗勞德數Fr變化的模擬結果與伯努利方程理論結果對比如圖5所示，可以看出二者十分吻合．

表1 模型計算條件Tab.1 Computational conditions of model

圖4 垂向平均流速模擬結果與實測結果對比Fig.4 Comparison of mean vertical velocity between numerical simulation and measurement results

圖5 樁前水位壅高值模擬結果與理論結果對比Fig.5Comparison of water level rising in front of the pile between numerical simulation and theoretical results

3 斷面形狀對繞流流場的影響

建立20,000,m×18.6,m的矩形數值水槽，計算區(qū)域左端為來流邊界，右端為出流邊界，上、下側均為滑移固邊界．底坡J＝1/100,000，槽底糙率0.015，水深20,m．均勻流流速為u0＝1.55,m/s．分別模擬斷面形狀為圓形(D＝0.6,m)、正方形(邊長L＝0.6,m)和正菱形(對角線長B＝0.6,m)的樁的繞流流態(tài)．樁置于水槽中心(見圖6)，計算區(qū)域及網格劃分見圖7.

圖6 不同斷面形狀樁繞流平面布置示意Fig.6Sketch of plane layout of flow around piles with different cross sections

不同斷面形狀樁的流速場模擬結果對比如圖8所示，水位模擬結果對比如圖9所示，橫斷面與縱剖面流速u與平均流流速u0之比u/u0的計算結果對比如圖10所示，橫斷面與縱剖面水位變化相對值(HH0)/D計算結果對比如圖11所示．從結果可以看出，在橫向上，圓樁與正方形樁及正菱形樁的流速分布不同，貼近樁體處圓樁流速最大，但隨后有所下降，在遠離樁體處圓樁的流速最小．正方形樁與正菱形樁相比，由于斷面形狀略大，對水流的擠壓作用也略大．從縱向看，正菱形樁后流速最小，圓樁與正方形樁后的流速較接近，表明其阻流效果好于圓樁和正方形樁．原因與正菱形樁的形狀有關．

圖7 不同斷面形狀的樁繞流計算區(qū)域及網格劃分(單位：m)Fig.7 Computational domain and mesh discretization of flow around piles with different cross setions(unit：m)

圖8 不同斷面形狀樁的流速場模擬結果對比(單位：m/s) Fig.8 Comparison of simulation results of flow velocity field for piles with different cross sections(unit：m/s)

圖9 不同斷面形狀樁的水位模擬結果對比(單位：m)Fig.9 Comparison of simulation results of water surface elevation for piles with different cross sections(unit：m)

圖10 不同斷面形狀的樁繞流u/u0計算結果對比Fig.10 Comparison of u/u0of flow around piles with different cross sections

圖11 不同斷面形狀的樁繞流(H-H0)/D計算結果對比Fig.11 Comparison of (H-H0)/D of flow around piles with different cross sections

4 透空式阻流堤阻流效果討論

工程中直接采用不透水堤來調整流速和流態(tài)往往會引起局部水流的劇烈改變，給工程帶來各種不利影響，如由于流速劇烈改變引起泥沙回淤或局部流速過大造成船舶航行或靠泊困難等．而采用透空式群樁阻流堤既可以達到調整流速的目的，又不至于引起當地流場的劇烈改變，作用效果比較和緩，但為了確定合理的群樁布置方案，需要了解群樁阻流效果．以下針對不同排列形式的樁排，研究其阻流效果．

選擇單排并列、兩排并列和3排并列3種布置形式進行研究，其中堤長按10根樁并列考慮．

為了嚴格控制水位和流速條件，數值水槽的上下游開邊界均給恒定水位．流速初值賦0．經試算，達到指定流速值需要的水槽長度應大于500,km，故實際取500,km，水槽寬100.8,m．水槽兩側邊界設置為自由滑移邊界．底坡J＝1/100,000，水深20,m，槽底糙率0.015，均勻流流速u0＝1.55,m/s．樁橫斷面選擇直徑D＝0.6,m的圓樁和邊長L＝0.6,m的方樁進行對比，橫向、縱向樁心距SH＝SV變化范圍為(2～4)D(L).為了減少網格數量，樁的上下游兩側大部分區(qū)域均采用矩形網格，樁周圍采用三角形網格，之間平滑過渡(見圖12)．多排樁按等間距梅花排列，見圖13．

圖12 圓樁樁排計算區(qū)域及網格劃分示意(單位：m)Fig.12Computational domain and mesh discretization of cylindrical pile rows(unit：m)

圖13 梅花排列示意Fig.13 Sketch of staggered array of pile group

由計算結果統(tǒng)計各種布置方式樁后(30～200)· D(L)處流速折減率處為上限，200,D(L)處為下限)、堤長及百米樁數如表2所示．圖14和圖15分別給出了圓樁和方樁在不同樁心距和排列情況下樁后不同位置處流速相對值的橫向分布結果，圖中y為沿樁排軸線方向坐標，y/D為相對距離．可以看出，方樁的阻流效果略好于圓樁.

表2 不同樁心距和排列群樁阻流堤后阻流效果數值模擬結果Tab.2Simulation results of the resistance effect in the rear side of pile-group with different spacings and array layouts

圖14 不同樁間距和排列圓樁的阻流效果Fig.14 Resistance effect of cylindrical pile group with different spacings and array layouts

圖15 不同樁心距和排列方樁的阻流效果Fig.15 Resistance effect of square pile group with different spacings and array layouts

5 群樁概化模擬

某碼頭為高樁墩式結構，碼頭軸線與流場中最大流速方向呈17°角，碼頭由17個墩臺構成，考慮到計算工作量，本研究暫模擬其東端4個墩臺．墩臺尺寸為31,m×26,m，由39根直徑D＝1.4,m、不規(guī)則排列的圓樁(見圖16)支撐．

矩形數值水槽長40,000,m、寬1,000,m．群樁位于水槽中央，水槽上下兩側為固邊界，與群樁有足夠的間距，水槽左、右端為開邊界．水槽底坡J＝1/100,000，水深25,m，槽底糙率0.015，均勻流流速u0＝1.80,m/s．

圖16 碼頭群樁平面布置及群樁周圍細網格剖分Fig.16Plane layout of pile group of wharf and fine mesh discretization

模擬分為直接和概化兩種情況．直接模擬是將單樁外輪廓當作固邊界對計算區(qū)域進行剖分(見圖16)，而概化模擬是將整個墩臺范圍劃分為8個單元

(見圖17)，群樁不等距分布在8個單元中(見圖18). 直接模擬單元總數為24,564，計算時間步長為0.01,s，概化模擬單元總數為1,852，計算時間步長為0.3,s．概化模擬計算量相當于直接模擬計算量的近1/400，計算效率大大提高．

圖17 群樁概化模擬粗網格剖分Fig.17 Coarse mesh for simplification simulation of pile group

圖18 概化模型中群樁分布Fig.18 Real pile group array of simplification model

從解決工程問題的角度出發(fā)，概化模擬中把不規(guī)則排列的39根樁簡化成40根8×5矩形排列的群樁，橫向樁心距SH＝3D，縱向樁心距SV＝4D．根據筆者利用計算流體力學軟件FLUNET對不同樁數和樁間距的矩形排列群樁阻力系數進行數值模擬計算的結果[16]，本算例中繞流阻力校正系數取40根樁分布在8個單元中，每個單元中有5根樁，群樁所在單元繞流阻力校正系數．根據樁的雷諾數，單樁阻力系數取0.4．

群樁所在單元面積為101,m2，阻流面積為95,m2，由式(11)算得折合水深ΔHt=23.5 m．

樁群概化模擬與直接模擬流場分布對比如圖19所示，圖中Δh為群樁引起的水位變化，負號表示水位降低．可以看出，二者流速是相符的，水位誤差在2%～5%．結果表明，這種考慮群樁阻力和遮流效果的概化方法既可有效提高計算效率，又有較好的計算精度．

圖19 群樁概化與直接模擬流場分布對比Fig.19Comparison of the flow field of simplification simulation for pile group against direct simulation

6 結 論

(1) 采用COASTALTOOL海岸軟件中的潮流模塊TC2D模擬了樁繞流問題，結果表明樁前水位壅高計算結果與理論值吻合較好，流速和水槽試驗實測結果一致，樁后較小范圍內的復雜流場模擬結果略差，但尾流區(qū)以外計算結果與試驗結果是吻合的．

(2) 通過對比不同斷面形狀的樁繞流模擬結果發(fā)現，在橫向上圓樁對水流的擠壓效果最明顯，而在縱向上正菱形樁的阻流效果好于圓樁和正方形樁，這應與樁的斷面形狀有關．

(3) 提出了考慮底部摩阻力及樁群繞流阻力的綜合阻力計算公式，公式考慮了群樁中各樁間的遮蔽作用；提出了初步考慮樁間相互影響的群樁概化方法. 群樁概化與直接模擬結果的對比表明，本文提出的群樁計算概化模式具有較好的精度．

(4) 研究了透空式群樁阻流堤斷面形狀、排列方式以及樁心距等因素對阻流效果的影響，得到了不同排列群樁后的流速折減率，結果表明正方形樁的阻流效果略好于圓樁．

［1］ Zdravkovich M M. Review of flow interference between two circular cylinders in various arrangements[J]. Journal of Fluids Engineering，1977，99(4)：618-633.

［2］ Williamson C H K. Evolution of a single wake behind apair of bluff bodies[J]. Journal of Fluid Mechanics，1985，159(1)：l-18.

［3］ 李景銀，Lam K，Chan K T，等. 繞正方形排列的順排的四個圓柱的流動研究[J]. 工程熱物理學報，2004，25(1)：59-62.

Li Jingyin，Lam K，Chan K T，et al. Stydy on the cross flow around four cylinders in an in-line square arrangement at low Reynolds numbers[J]. Journal of Engineering Thermophysics，2004，25(1)：59-62(in Chinese).

［4］ 鄧紹云，張嘉利. 樁群阻力測試的研究[J]. 華北水利水電學院學報，2007，28(2)：86-90.

Deng Shaoyun，Zhang Jiali. Study on experiment of testing for drag force of water flow around groups of piles[J]. Journal of North China Institute of Water Conservancy and Hydroelectric Power，2007，28(2)：86-90(in Chinese).

［5］ 鄧紹云. 樁柱水流繞流阻力特性及其計算[J]. 中國港灣建設，2007(1)：14-21.

Deng Shaoyun. Drag force characteristics and calculation of water flow around pile[J]. China Harbour Engineering，2007(1)：14-21(in Chinese).

［6］ 黃筱云，李紹武. 自適應四叉樹網格下的N-S方程數值求解模型[J]. 天津大學學報：自然科學與工程技術版，2013，46(1)：58-66.

Huang Xiaoyun，Li Shaowu. Numerical N-S equation solver based on adaptive ouadtree mesh[J]. Journal of Tianjin University：Science and Technology，2013，46(1)：58-66(in Chinese).

［7］ Stansby P K. A numerical study of vortex shedding from one and two circular cylinders[J]. Aeronautical Quarterly，1981，32(1)：48-71.

［8］ 吳文權. 流體繞多個鈍體不穩(wěn)定分離流動數值仿真[J]. 華東工業(yè)大學學報，1997，19(3)：1-8.

Wu Wenquan. Numerical simulation for unstable separated flow past multi-blunts[J]. Journal of Huadong Poly-technic University，1997，19(3)：1-8(in Chinese).

［9］ Farrant T，Tan M，Price W G. A cell boundary element method applied to laminar vortex-shedding from arrays of cylinder in various arrangements[J]. Journal of Fluids and Structures，2000，14(3)：375-402.

［10］ 張愛社，張 陵. 等邊布置三圓樁繞流的數值分析[J]. 應用力學學報，2003，20(l)：142-150.

Zhang Aishe，Zhang Ling. Numerical simulation of three equispace dcircular cylinders[J]. Chinese Journal of Applied Mechanics，2003，20(l)：142-150(in Chinese).

［11］ 唐士芳. 二維潮流數值水槽的樁群數值模擬[J]. 中國港灣建設，2002，6(3)：14-21.

Tang Shifang. Numerical simulation for pile group innumerical water flume of two dimensional tidal flow[J]. China Harbour Engineering，2002，6(3)：14-21(in Chinese).

［12］ 唐士芳，李 蓓. 樁群阻力影響下的潮流數值模擬研究[J]. 中國港灣建設，2001，10(5)：25-29.

Tang Shifang，Li Pei. Study on numerical simulation of tidal flow influenced by pile group resistance[J]. China Harbour Engineering，2001，10(5)：25-29(in Chinese).

［13］ 解鳴曉，張 瑋，謝慧姣. 樁群數值模擬中的概化方法研究[J]. 水動力學研究與進展，2008，A23(4)：464-471.

Xie Mingxiao，Zhang Wei，Xie Huijiao. Simplification method in numerical modeling of bridge pier group[J]. Chinese Journal of Hydrodynamics，2008，A23(4)：464-471(in Chinese).

［14］ 李紹武，盧麗峰，時 鐘. 河口準三維涌潮數學模型研究[J]. 水動力學研究與進展，2004，A19(4)：407-415.

Li Shaowu，Lu Lifeng，Shi Zhong. A quasi-3D numerical model for estuarinetidal bore[J]. Chinese Journal of Hydrodynamics，2004，A19(4)：407-415(in Chinese).

［15］ Roulund A，Sumer B M，Freds?e J，et al. Numerical and experimental investigation of flow and scour around a circular pile[J]. Journal of Fluid Mechanics，2005，534：351-401.

［16］ 莊 茜. 建筑物繞流數學模型理論及應用研究[D]. 天津：天津大學建筑工程學院，2013.

Zhuang Qian. Study on Theory and Application of Mathematical Model for Flow Around Structures[D]. Tianjin：School of Civil Engineering，Tianjin University，2013(in Chinese).

(責任編輯：樊素英)

Numerical Simulation for Flow Around Pile Group Based on FVM

Zhuang Qian1,2,3，Li Shaowu1，Qi Zepeng1
(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China；2. School of Civil Engineering，Tianjin Chengjian University，Tianjin 300384，China；3. Key Laboratory of Water and Sediment Science and Water Hazard Prevention，Changsha 410114，China)

Flow around pile is studied by using the tide module TC2D of COASTALTOOL software. The effects of pile transection shape，array layout and center distance of pile group on the flow around piles are investigated，respectively，by adopting fine grid model. As for the coarse grid model，a calculation formula is proposed for the resistance of pile group inside an element based on the superposition principle of resistance. A simplification method is put forward to account for the shielding effect of pile group in an element by introducing the equivalent water depth for the sheltered area of the piles，by which the shielding effect among piles of a pile group is considered.

flow around pile；resistance of pile group；shielding effect；simplification

TV131

A

0493-2137(2015)05-0445-10

10.11784/tdxbz201403091

2014-03-27；

2014-09-29.

國家自然科學基金資助項目(51379143)；國家創(chuàng)新研究群體基金資助項目(51021004)；水沙科學與水災害防治湖南省重點實驗室基金資助項目(2014SS01).

莊 茜（1980— ），女，博士研究生，講師，zq@tcu.edu.cn.

李紹武，lishaowu@tju.edu.cn.