孤立波作用下排柱附近流動特性數值分析

蔣昌波，鄧 涯，姚 宇，鄧 斌

（1.長沙理工大學 水利工程學院，湖南 長沙410004；2.水沙科學與水災害防治湖南省重點實驗室，湖南 長沙410004）

海嘯是一種破壞性極強的海洋災害.海嘯波首波達到近海岸時波高急劇增加，非線性特征增強，近似于孤立波，所以對海嘯通常采用孤立波來模擬［1］.在海嘯波向近岸推進過程中，水體內部蘊含著巨大的能量，會對近岸建筑物產生巨大沖擊和潛在破壞［2］，海嘯的危害在2004年印度洋海嘯和2011年日本海嘯中已經得到充分表現.因而對于海嘯波與海洋建筑物相互作用的研究已經成為海嘯災害防治相關研究的熱點.近年來隨著海洋開發的不斷拓深，多圓柱結構開始廣泛應用于各種海洋和海岸工程建筑物，如海底管線、樁式防波堤、海洋鉆井平臺和跨海大橋的基礎等.所以研究孤立波與多圓柱之間的作用機理，對于評估海嘯波對這類建筑物的破壞作用具有十分重要的指導意義.

由于圓柱間水流結構的相互影響，多圓柱繞流問題相比傳統的單圓柱繞流問題更為復雜.在孤立波與多圓柱相互作用過程中，自由液面變化劇烈，圓柱附近水體的非恒定特性和三維特征明顯.目前國內外對此類問題的研究較少且以實驗研究和數值模擬為主.Huang等［3］以排樁式防波堤為背景對孤立波與圓形排柱間的相互作用進行了物理模型實驗研究，分析了排柱間隙大小與孤立波波高衰減之間的規律.Liu［4］在Huang［3］單排圓柱實驗的基礎上，對孤立波波高的衰減和排柱附近的流場進行了數值模擬研究.但由于受到淺水波方程的限制，該模型不適應于非線性較強的孤立波，也不能描述沿水深方向的流動特征.

基于Navier-Storks方程的OpenFOAM 作為CFD 開源程序包，以易擴展、并行計算穩定高效和求解方法先進等優點在計算流體力學領域得到了廣泛的應用，近年來也開始應用于波浪與圓柱的相互作用問 題.Mo 等［5］和Lara 等［6］分 別 采 用Open-FOAM 開源程序包對孤立波與三圓柱間相互作用進行了數值模擬研究，分析了三圓柱附近的自由液面、流速、和動水壓強的變化規律.但由于研究側重點不同，Mo等［5］和Lara等［6］并沒有分析多圓柱附近渦結構的發展和流場特性.國內對孤立波與圓柱間相互作用的數值模擬研究甚少，主要為規則波作用下的單個圓柱模型研究，如祝兵［7］和曹洪建［8］.

綜上所述，為了彌補孤立波與排柱相互作用的研究中現有數值方法（如淺水波方程）的不足，分析排柱附近垂向上的水流結構，本文擬在Huang等［3］圓形排柱物理實驗研究的基礎上，利用OpenFOAM中的不可壓縮氣液兩相流求解器interFoam 來建立三維數值波浪水槽，模擬孤立波與單排圓柱間的相互作用過程.通過計算分析排柱附近自由液面、平面／立面流線以及平面渦量場的變化來研究排柱附近的三維流動特性.

1 數值模型

1.1 控制方程

數值模型是求解連續不可壓縮流體的氣液兩相流雷諾時均N-S 方程.自由面的捕捉是采用VOF方法，引入了流體體積分數φB，其值大小代表液體體積在計算網格中所占的比例，φB＝1 表示液體，φB＝0表示氣體，介于0和1之間表示該網格為氣液交界面，其輸運方程如下：

式中：u為速度，ur為相對速度.該方法與經典VOF模型中的體積分數輸運方程相比多了等式左邊的最后一項，稱為人工壓縮項，其作用是使其對自由液面的捕捉更加精確.

氣液兩相流的動量方程中在方程右端增加了包含有體積分數的表面張力一項［9］，相應控制方程如下：

連續性方程：

動量方程：

式中：u＝［u，v，w］是笛卡爾坐標系下流速場；p＊為動水壓強；ρ為密度；g 為重力加速度；μ 為流體的動力黏度；ρ＝ρ（x）大小與計算網格內氣液所占體積有關；τ為雷諾應力：

式 中：μt 為 湍流黏 度；k 為 紊 動 動 能；I為Kronecker張量；S 為應變速率（1／2（▽u＋（▽u）T））；▽為向量算子（?／?x，?／?y，?／?z）T；x＝（x，y，z）為Cartesian坐標系.式（3）右端最后一項σTκ ▽φB為表面張力，σT為表面張力系數，取值為0.074kg／s2；κ 為界面曲率，即

OpenFOAM 使用有限體積法對N-S方程進行離散，方程求解是采用PIMPLE 算法實現壓力和速度的解耦，時間離散采用Euler格式，壓力梯度離散采用Gauss線性離散格式，拉普拉斯項離散采用修正后Gauss線性離散格式，具體離散格式與過程詳見文獻［10］.

1.2 數值造波

1.2.1 造波和消波方法 采用Jacobsen［11］提出的修正速度入口方法進行造波.與通常的速度入口造波方法不同，該方法中造波邊界的類型分為干網格、濕網格和臨界網格3種，其中干網格的邊界條件為

式中：n為邊界網格面的單位法向向量.濕網格邊界條件（u和φB）是根據相應的波浪理論計算所得.模型對臨界網格進行了濕面修正.波浪理論計算的自由液面曲線與網格邊界存在2個交點，取兩點相連直線作為臨界網格的濕面邊界，得到網格濕面部分的中心位置，以此計算出臨界網格的入流速度和體積分數.對臨界網格的處理，避免了入口處自由液面的虛擬振蕩.

消波方法采用松弛因子法，對消波區網格的計算變量φ（u和φB）進行數值修正：

式中：αR為松弛因子，消波區域和非消波區域交界面上αR＝1，在消波區域αR∈（0，1），φcal為變量的計算值，φtar為變量的消減目標值.

1.2.2 造波驗證 本文對數值波浪水槽進行了孤立波造波和傳播的驗證.驗證水槽長11 m、高0.3m，寬為一個網格大小，兩端設置長1m 的數值消波區（Relexation Zone）.模擬水深0.15 m，孤立波波高0.07m，造波起點設在x＝3m 位置（x＝0為計算左端邊界），計算網格大小為0.01m，時間步長為自適應調整步長，驗證結果如圖1所示.孤立波波面和流速理論公式為

圖1 孤立波造波驗證Fig.1 Verification of solitary wave generation and propagation

式 中：η 為 自 由 液 面 變 化 值，H 為 波 高，h 為 水 深，c為波速.

如圖1（a）所示數值生成的孤立波波形與理論波形一致，其中η＊為自由液面高程值，數值計算的波面歷時數據取自數值造波點后3m 處（x＝6m）.孤立波在傳播過程中，除造波點后2m 區域內（x＝3～5m）波高有略微振蕩外，在其他區域波高基本保持不變，水槽兩端的消波效果良好（圖1（b）中灰實線所示）.通過對數值模擬的相鄰時刻孤立波水平傳播距離進行分析，得到孤立波的傳播速度ccal＝1.43m／s，根據式（8）計算理論波速為ctho＝1.47 m／s，兩者較為符合（誤差小于3%）.因而該模型可以有效的模擬孤立波的生成和傳播過程.

2 模型設置與驗證

為了驗證數值模型能否合理地模擬孤立波與排柱相互作用的過程，數值模型的設置與Huang等［2］的物理模型實驗一致，其模型設置如圖2（a）所示，以孤立波傳播方向為x 方向，水槽斷面方向為y 方向，水深方向為z方向建立空間直角坐標系，原點位于入口端水槽邊壁的底部位置.水槽寬0.54 m，排柱模型直徑d＝0.03 m，相鄰圓柱圓心間距S＝0.045m，在x＝8m 處沿水槽橫斷面方向對稱布置12根，在排柱模型前3m 和后1.5m 位置分別布置有浪高儀G1、G2.

圖2 數值水槽和排柱附近網格劃分Fig.2 Wave numerical tank and computational grids around cylinders

為節省計算資源，數值計算域選定水槽中間寬0.27m 的區域，包含有6根圓柱.上下選用對稱邊界以消除邊壁效應的影響，最外層圓柱圓心距邊界距離為S／2，如圖2（a）陰影區域所示.數值水槽在長度和高度方向的設置與圖1（b）所示相似，造波起點位于排柱前5m 位置，水槽入口端和出口端設置有長為1m 的數值消波區.同時為了保證排柱前后入射波和透射波的穩定傳播，排柱前后分別設置有5和2m的富余長度，總計算區域長為11m.數值水槽的計算網格采用分區多塊結構網格.計算區域中y 和z 方向網格大小相等dy＝dz＝0.006m，x 方向網格大小則是分區設置.在排柱前1m 至后0.5m區域內，x 方向網格大小為dx＝0.006m，在造波起點和排柱后2m 位置點dx＝0.01m，在入口和出口位置點dx＝0.02m，相鄰位置點間網格進行均勻過渡.在水槽網格的基礎上，排柱附近網格采用snappyHexMesh工具［12］進行劃分：對圓柱表面采用單層貼體網格，排柱附近由加密的結構網格填充，與外層結構網格的銜接則是采用三角棱柱非結構網格（如圖2（b）所示），網格總數為2 908 350.進一步加密圓柱附近的網格后計算結果改進不明顯，說明排柱附近的網格設置可以滿足計算要求.同時本文進行了單柱模型的對比模擬，單柱附近網格生成方法與排柱情況相同，網格總數為2 839 600.

本文選用一個典型的實驗工況（水深h＝0.15 m，孤立波波高H＝0.07m）進行研究，湍流模型采用文獻［11］建議的SSTk-ω 模型，時間步長為自適應時間步長，計算時長為8s，以保證反射波傳至G1測點.

本模型中具體的邊界條件設置為：水槽兩側為對稱邊界，底部和圓柱表面為固壁無滑移邊界，上方空氣為自由進出流邊界，其中空氣區域高0.15 m，以防止水流溢出.

圖3 自由液面歷時曲線的數值計算和實驗結果對比Fig.3 Comparison between numerical and experimental time series of free surface elevation0.86，Kd＝0.50，實驗和數值結果符合較好.

3 結果分析與討論

3.1 自由液面的變化

圖4 不同時刻的自由液面Fig.4 Free surfaces at different time steps

圖5 各個圓柱附近的自由液面變化Fig.5 Surface elevations around cylinders

與傳統的沿水深積分的二維數值波浪模型（如淺水波方程）相比，該模型能更合理的模擬排柱附自由液面的三維特性：圖5表明單個圓柱表面自由液面壅高沿中軸線（θ＝π）兩側近似對稱分布.在柱前位置（θ＝π）壅高達到最大，自由液面變化相對值為1.7，在柱后附近位置（θ＝π／8（15π／8））達到最小，自由液面變化相對值為0.6，兩者差值的相對值為1.1.在θ＝π到θ＝π／2（3π／2）位置，壅高逐漸減小，液面變化相對值為0.1.在θ＝π／2（3π／2）到θ＝π／4（7π／4）位置，壅高出現了急劇下降，液面變化相對值0.6.在θ＜π／4，θ＞7π／4范圍內，壅高變化幅度較小.

3.2 平面和立面流線特征

孤立波與排柱相互作用過程中，水流結構復雜多變，三維特征明顯.故本文分別選取z／h＝0.83位置x-y 平面（平面的選取保證其始終位于自由液面以下，并且盡可能接近靜水面位置z／h＝1）和y／d＝9位置x-z 立面（立面的選取位于圓柱3和圓柱4的間隙中心線上）為例，來研究排柱附近水流的三維流線特征，如圖6 所示為不同時刻排柱附近的流線.圖6中（a）～（e）為不同時刻的平面流線圖，（f）～（i）為不同時刻的立面流線圖，平面和立面位置如圖中點畫線所示.

3.3 水平渦量場分析

排柱結構除了利用樁柱自身的反射外，還依靠排柱間隙產生的高速水流和柱后的尾渦脫落來增強入射波波能的耗散.為研究排柱附近自由液面渦動強度的變化規律，本文選取z／h＝0.83 位置的x-y平面（與3.2節一致）作為研究對象.由于各個圓柱附近水流結構相似，為了清晰展示排柱附近的流場結構和渦量特征，文中僅選取中間3根圓柱附近的區域作為研究對象，并且與相同孤立波作用下單柱的情況進行了比較，如圖7所示.

圖6 平面（z／h＝0.83）和立面（y／d＝9）不同時刻流線Fig.6 Streamlines on the horizontal（z／h＝0.83）and the vertical（y／d ＝9）planes of at different time steps

圖7 z／h＝0.83位置平面上渦量（s-1）Fig.7 Vorticity on plane of z／h＝0.83（s-1）

4 結 語

本文利用OpenFOAM 開源程序包建立三維數值波浪水槽，研究了孤立波與排柱相互作用下排柱附近流動特性的演變規律.圓柱表面邊界層網格與外層結構網格采用了三角棱柱網格進行銜接.與傳統的二維數值波浪模型相比，該模型能更合理的模擬排柱附近波浪傳播的三維特性.研究結果表明：孤立波穿過排柱時，自由液面顯著壅高，并沿柱面呈三維分布，隨后在排柱前后分別形成沿上游傳播的反射波和沿下游傳播的透射波.排柱中各圓柱附近流動特性相似，柱前壅高最大時柱間水體以水平運動為主，垂向流速較小，具有典型的淺水長波流動特征；圓柱下游形成對稱分布的漩渦，并逐漸耗散消失.與孤立波作用下單柱附近的渦量特征相比，排柱后方渦動強度較大，且尾渦的發展縱向相對伸長、橫向相對收縮.本文的研究成果將為進一步分析排柱的受力特征以及排柱附近的泥沙運動規律提供參考.

（

）：

［1］CHANG Y H，HWANG K S，HWUNG H H.Largescale laboratory measurements of solitary wave inundation on a 1：20slope［J］.Coastal Engineering，2009，56（10）：1022-1034.

［2］陳杰，蔣昌波，鄧斌，等.海嘯作用下岸灘演變與床沙組成變化研究綜述［J］.水科學進展，2013：750-758.CHEN Jie，JIANG Chang-bo，DENG Bin，et al.Review of beach profile changes and sorting of sand grains by tsunami waves ［J］.Advances in Water Science，2013，24（5）：750-758.

［3］HUANG Z，YUAN Z.Transmission of solitary waves through slotted barriers：A laboratory study with analysis by a long wave approximation［J］.Journal of Hydro-Environment Research，2010，3（4）：179-185.

［4］LIU H，GHIDAOUI M S，HUANG Z，et al.Numerical investigation of the interactions between solitary waves and pile breakwaters using BGK-based methods［J］.Computers ＆ Mathematics with Applications，2011，61（12）：3668-3677.

［5］MO W，LIU P L F.Three dimensional numerical simulations for non-breaking solitary wave interacting with a group of slender vertical cylinders［J］.International Journal of Naval Architecture and Ocean Engineering，2009，1（1）：20-28.

［6］LARA J L，HIGUERA P，GUANCHE R，et al.Wave interaction with piled structures：application with IHFOAM［C］∥ASME 2013 32nd International Conference on Ocean，Offshore and Arctic Engineering.Nantes，France：American Society of Mechanical Engineers，2013：V007T08A078.

［7］祝兵，宋隨弟，譚長建.三維波浪作用下大直徑圓柱繞流的數值模擬［J］.西南交通大學學報，2012，47（2）：224-229.ZHU Bing，SONG Sui-di，TAN Chang-jian.Numerical simulation for diffraction around large-diameter circular cylinder subjected to three-dimension wave［J］.Journal of Southwest Jiao Tong University，2012，47（2）：224-229.

［8］CAO H J，ZHA J J，WAN D.Numerical simulation of wave run-up around a vertical cylinder［C］∥Proceedings of the 21st International Offshore and Polar Engineering Conference.Maui，Hawaii，USA：International Society of Offshore and Polar Engineers，2011：726.

［9］RUSCHE H.Computational fluid dynamics of dispersed two-phase flows at high phase fractions［D］.London，UK：University of London，2002.

［10］Open CFD.OpenFOAM User Guide，Version 2.2.1［EB／OL］.2013-10-20.http：∥www.openfoam.com／.

［11］JACOBSEN N G，FUHRMAN D R，Freds?e J.A wave generation toolbox for the open‐source CFD library：OpenFoam［J］.International Journal for Numerical Methods in Fluids，2012，70（9）：1073-1088.

［12］JACKSON A.A comprehensive tour of snappyhexmesh［R］.Darmstadt，Germany：The 7th Open-FOAM Workshop，2012.