透水魚礁型潛堤內(nèi)部流場及繞流特性模擬

匡翠萍， 鄭宇華， 顧杰， 韓雪健

（1.同濟大學 土木工程學院，上海 200092；2.上海海洋大學 海洋生態(tài)與環(huán)境學院，上海 201306）

近年來，在國內(nèi)外大力提倡效法自然的工程設計理念背景下，新興的環(huán)境友好型海岸防護建筑物——魚礁型潛堤應運而生。魚礁型潛堤將人工魚礁的設計融合到離岸潛堤的結構中，這一概念最早由Armono和Hall［1］提出并通過試驗證明了這兩類構筑物結合設計使用的可行性。我國一直致力于提升海域、海島和海岸帶的環(huán)境和生態(tài)價值。自2010年起，我國京津冀旅游勝地北戴河為修復侵蝕嚴重的海岸開展了一系列的災害治理工程，工程岸線長度3.5km，工程后海岸修復成效顯著。值得注意的是，在工程中首次采用了魚礁型潛堤代替了傳統(tǒng)的堆石離岸堤應用于海岸恢復治理，改善了近岸海域的生態(tài)環(huán)境，增加了漁業(yè)資源［2］。魚礁型潛堤多設計為空心的透水結構，投放于工程海域中礁體內(nèi)部和周圍會形成連續(xù)多變的流場形態(tài)，能夠迎合不同海洋生物對水流條件的趨向性［3］，兼具生態(tài)修復和改造流場的效果。魚礁型潛堤形成的流場本質(zhì)上體現(xiàn)了復雜的鈍體繞流現(xiàn)象，包含了高度非線性的湍流流動過程，如流體分離和渦旋演化，這類問題往往難以得到解析解，卻又是準確評估魚礁型潛堤產(chǎn)生的流場效應和生態(tài)效應的重要依據(jù)［4］。

至今，已有眾多學者利用物模試驗和數(shù)值模擬方法研究了魚礁型潛堤繞流流場特性。吳建［5］基于ADV測流結果分析了不同來流速度下方體礁型潛堤周圍的流場變化，分別采用重整化群k-ε（Renormalization Groupk-ε）和可實現(xiàn)k-ε（Realizablek-ε）湍流模型模擬計算了單體箱型人工魚礁的繞流特征，結果表明2個湍流模型對礁體周圍的渦旋結構都有一定的預測能力，而可實現(xiàn)kε模型在流線曲率變化較大的地方比重整化群k-ε更準確一些。Li等［6］通過重構的k-ε湍流數(shù)學模型和PIV測流技術分析了圓管型人工魚礁周圍上升流的強度和尺度，發(fā)現(xiàn)上升流的強度和尺度隨著魚礁數(shù)量的增多而增大，當礁高與水深比值為0.224時，魚礁組合產(chǎn)生的效應更強。王者也和李爽［7］使用并行大渦模擬模式（Parallelizeda Large-eddy Simulation Model， PALM）及示蹤物模型，研究了不同流速條件下方型魚礁型潛堤對流場形態(tài)和湍流動能收支的影響。崔勇等［8］和胡聰?shù)龋?］均采用水槽試驗和重整化群k-ε湍流模型分析了不同縱橫布設間距對圓筒形礁體流場的影響，得到上升流體積隨礁體橫縱向布設變化的規(guī)律。Tang等［10］根據(jù)氫氣球試驗結果對比了3個雙方程k-ε模型（標準k-ε、重整化群k-ε和可實現(xiàn)k-ε模型）對方體人工魚礁流場的適用性，發(fā)現(xiàn)標準k-ε（Standardk-ε）模擬結果與試驗得到的礁體外流場中的流線變化基本一致。Maslov等［11］以不規(guī)則多功能人工魚礁為研究對象，將雙方程k-ε、k-ω（ω為比耗散率［12］，μt=ρk/ω）模型以及剪切應力輸運模型（Shear Stress Transport， SST）模擬得到的礁體周圍流速分布與ADV測量結果進行對比，發(fā)現(xiàn)以單調(diào)迎風格式離散的k-ω模型計算結果與實測值吻合最好。Kim等［13］基于粒子圖像測速（particle image velocimetry，PIV）的試驗結果，比較了標準k-ε、重整化群k-ε、k-ω、SSTk-ω和雷諾應力（Reynolds Stress Model，RSM）模型在半球形人工魚礁流場模擬中，礁體尾流的規(guī)模和位置的不同，發(fā)現(xiàn)RSM模型能夠最準確地預測魚礁尾流的規(guī)模和位置，然而RSM模型的計算時間明顯高于其他另外4種湍流模型，計算量大且難以收斂，因此在工程中不常采用。近年來越來越多學者關注魚礁型潛堤的繞流流場特性，大多數(shù)仍通過改變外界條件研究礁體周圍水流結構（如上升流）的變化，而對于透水魚礁型潛堤流場中湍流擾動最復雜的部分，即礁體內(nèi)部的流動分離特征尚未有詳細討論。針對魚礁型潛堤的數(shù)值計算，由于每種湍流模型基于的各種假設和包含的經(jīng)驗常數(shù)都具有一定局限性，可以看到各家研究中湍流模型的采用不單一且不唯一，根據(jù)特定水流條件下結構物的不同，湍流模型的選擇也可能不同。從已有的研究成果可以看出，雙方程湍流模型已經(jīng)能夠基本刻畫魚礁型潛堤周圍的紊流情況，其中以雙方程標準k-ε和重整化群k-ε模型應用最為廣泛。隨著計算機技術的飛速發(fā)展和計算流體力學的精度要求不斷提高，大渦模擬也被更多地應用于結構-流體相互作用形成繞流流場的精細化研究，其在復雜的魚礁型潛堤繞流流場內(nèi)精確求解某個尺度上所有的湍流運動體現(xiàn)了很大的優(yōu)勢。可見，目前對于魚礁型潛堤周圍繞流流場的預測可選的模擬方式較多，可根據(jù)不同的環(huán)境條件選擇不同的湍流模型，而對于魚礁型潛堤內(nèi)部紊流場的刻畫方式目前仍未有足夠的研究成果支持統(tǒng)一定論。

本文以北戴河生態(tài)修復工程中投放的透水魚礁型潛堤為原型，通過合理比尺制作試驗模型，采用PIV技術對單體魚礁型潛堤模型進行測流試驗，分析魚礁型潛堤內(nèi)部和周圍形成的繞流流場特征?；谟邢摅w積法，建立與PIV試驗相同的數(shù)學模型，分別采用3個雙方程k-ε模型（標準k-ε、重整化群k-ε、可實現(xiàn)k-ε）和大渦模擬對魚礁型潛堤的繞流流場進行數(shù)值模擬，通過對比計算與實測得到的流場分布和礁體內(nèi)外的渦旋結構分析這4種湍流模型的準確性和適用性，為進一步了解和優(yōu)化其在海洋生態(tài)修復工程領域的應用提供參考。

1 水槽試驗

試驗在上海海洋大學水動力學實驗室的循環(huán)水槽中進行， 水槽試驗段長6.00m、寬0.45m、高0.55m，水槽底板及側壁均為鋼化玻璃，便于開展PIV流速測量，具體如圖1a所示。水槽系統(tǒng)采用循環(huán)方式供水，進口處底部為蜂窩狀鋼板，具有整流和消除大尺度渦結構的作用，進水段壁面向直道段沿曲線收縮，使得水流流速沿下游方向單調(diào)增加，保證了水體進入試驗段后流態(tài)平穩(wěn)。基于目標投放海域水深（約7m）和實驗水槽尺寸，將魚礁型潛堤按幾何比尺1：12.5進行模型設計，材料為亞克力有機玻璃，透光率為92%，模型長L=0.214m、寬B=0.152m、高H=0.224m，模型的頂板和迎流面上均開有4 個直徑為0.024m的圓孔，如圖1b所示。考慮到魚礁型潛堤的實際投放情況，北戴河海域平均流速約為0.3m·s-1，工程區(qū)內(nèi)流速基本小于0.10 m·s-1［14］，因此，在本試驗中，經(jīng)均勻流測試后，將進口流量設置為40 m3·h-1，相應入口平均流速為u0=0.091 4 m·s-1。魚礁型潛堤模型放置于水槽中央，距離入水口3m，水深為0.27m。試驗中雷諾數(shù)為18 121，屬高雷諾數(shù)的湍流。

圖1 PIV試驗及魚礁型潛堤模型示意Fig.1 Schematic diagram of PIV test stand and reef-type breakwater model

流速測量采用美國TSI公司研發(fā)的粒子圖像測速系統(tǒng)，連續(xù)記錄礁型潛堤近區(qū)瞬時速度場的時間序列數(shù)據(jù)，并計算出反映流場特性的速度剖面及相應的紊流脈動數(shù)據(jù)。PIV系統(tǒng)主要包括示蹤粒子、光學照明系統(tǒng)、圖像記錄系統(tǒng)和圖像處理系統(tǒng)4個部分。示蹤粒子采用聚氯乙烯微粉（PVC-6500），具有跟隨性好、化學性質(zhì)穩(wěn)定的優(yōu)點。照明光源采用雙脈沖激光器，片光厚約2 mm；圖像記錄系統(tǒng)主要通過高頻CCD相機（分辨率為2 048×2 048像素）連續(xù)曝光，在不同時段拍攝魚礁型潛堤近區(qū)流場的粒子圖像，采樣頻率為4.83Hz，即每秒獲得4.83對PIV圖像，進而得到4.83個速度場的數(shù)據(jù)。所得流場圖的后處理通過PIV系統(tǒng)自帶的Insight 4G軟件，基于互相關算法［15］對粒子速度場圖中的流速矢量進行擬合計算。研究表明，人工魚礁周圍的三維渦結構在由展向渦和流向渦組合而成，在垂向和水平向上差異顯著［16］。因此，試驗選取了3個不同位置的剖面進行流速測量，激光光源首先從水槽上方向下投射，激光垂直穿過礁體頂板的一列圓孔；隨后調(diào)整激光器位置，分別使激光水平穿過礁體迎流面的上、下2排圓孔，利用CCD相機依次捕捉激光照亮的粒子圖像得到3個剖面的流場圖，即圖1b魚礁型潛堤中的剖面（1）、（2）和（3）。每個剖面的測流每次連續(xù)捕捉500個瞬時流速場，然后將瞬時流場數(shù)據(jù)進行平均，得到時均流場圖，并獲得相應的流場參數(shù)。

2 數(shù)學模型

2.1 模型設置

模擬采用的魚礁型潛堤規(guī)格與PIV試驗中的一致。計算域尺寸的設定根據(jù)PIV試驗中水槽的尺寸和礁體潛堤模型的尺寸綜合決定。在數(shù)學模型中，將計算域的寬度設定為水槽寬度0.45m，計算域的高度設定為水深高度0.27m，計算域的長度則根據(jù)單體魚礁型潛堤的尺寸進行設定，為16倍的礁體長度（礁體前方5倍，礁體后方10倍），以確保礁體后方尾流區(qū)的范圍能達到10倍的礁體尺寸。

因多孔透水魚礁型潛堤結構較為復雜，采用非結構網(wǎng)格對魚礁型潛堤進行網(wǎng)格劃分，如圖2所示。魚礁型潛堤壁面和底部壁面與水流的交界處存在邊界層，邊界層很薄，但是其中的各項流體參數(shù)卻發(fā)生著劇烈的變化，存在較大的速度梯度。為準確模擬礁體近區(qū)的繞流流動，將近壁邊界層第1層網(wǎng)格厚度設為0.5mm，變化率為1.1，對魚礁型潛堤迎流面、礁頂及透水圓孔內(nèi)壁進行局部加密。數(shù)值計算結果的準確性很大程度上取決于網(wǎng)格的數(shù)量和質(zhì)量。為減小數(shù)值模擬過程中由于網(wǎng)格尺度所引起的誤差，首先針對網(wǎng)格數(shù)量無關性進行了研究，并保證網(wǎng)格質(zhì)量都在同一水平，通過改變模型表面和周圍網(wǎng)格的大小來確定最優(yōu)的網(wǎng)格數(shù)?？紤]到不同湍流模型對網(wǎng)格劃分及密度的要求不一樣，大渦模擬模型對網(wǎng)格要求最高，因此以大渦模擬模型計算結果為參考，依據(jù)Biron等［17］的方法進行網(wǎng)格無關性驗證，確定最合適的網(wǎng)格全局尺寸，即當不同網(wǎng)格模型的計算結果相近時，90%流場點的誤差控制在5%以內(nèi)，可以認為計算值與真實值比較接近。以最密網(wǎng)格（網(wǎng)格數(shù)量360萬）為基準計算不同網(wǎng)格數(shù)量下的標準差，結果如表1所示，網(wǎng)格在97萬時即已經(jīng)達到了網(wǎng)格獨立性要求，過密的網(wǎng)格對流場最高流速的精確度沒有明顯的提高?；诹鞴探唤缑娴牧鲃犹匦苑治觯C合考慮模型的計算時間和計算精度，采用網(wǎng)格總數(shù)為112萬時確定的全局網(wǎng)格尺寸。

表1 網(wǎng)格無關性驗證Tab.1 Verification of grid independence

圖2 數(shù)值水槽及魚礁型潛堤模型網(wǎng)格劃分示意Fig.2 Schematic diagram of computational domain and the grids distributed near a single reeftype breakwater

網(wǎng)格無關性標準差為

式中：Xi為第i次最大流速Umax；δ為第i次之后所有Xi的標準差；n為次數(shù)；g為重力加速度。

計算域邊界條件設置為：入口邊界條件選擇速度入口，流速根據(jù)水槽試驗的實際測量值給定； 出口邊界條件選擇壓力出口，平均靜壓P為零；魚礁型潛堤個體及水槽壁面均設置為無滑移固壁邊界，并采用無滑移邊界條件；由于魚礁型潛堤在淹沒狀態(tài)下水面較平穩(wěn)，自由水面采用剛蓋假定處理，設置為對稱邊界（在進行魚礁型潛堤結構周圍和內(nèi)部水力特性研究時，主要觀察其內(nèi)部和周圍的流場，而其水面起伏及水位變化對研究結果影響不大，考慮到計算資源及計算時間，采用剛蓋假定對其表面進行處理）。

根據(jù)魚礁型潛堤近區(qū)水流運動的特性，將流體假設為不可壓縮、定常、黏性流體，采用有限體積法離散控制方程，基于壓力求解器，壓力-速度耦合方式采用SIMPLE算法，壓力項采用二階方法離散，單元中心梯度插值方法采用基于單元體的最小二乘法插值。

2.2 湍流模型

魚礁型潛堤近區(qū)形成的繞流是空間上不規(guī)則和時間上無秩序的非線性湍流流體運動。湍流數(shù)值模擬可分為直接數(shù)值模擬和間接數(shù)值模擬，前者直接求解Navior-Stokes（N-S）方程得到有效的湍流信息，但由于N-S方程本身求解的復雜性和巨大的計算量，現(xiàn)有的計算資源往往難以滿足高雷諾數(shù)流動直接模擬的條件，因此限制了其應用范圍。為簡化計算過程使工程計算具有可行性，通常對湍流做某種程度的近似處理以簡化N-S方程進行間接數(shù)值模擬，主要手段包括雷諾平均法（Reynolds Average Navier-Stokes，RANS）和大渦模擬法（Large Eddy Simulation，LES）。RANS模型通過引入雷諾應力和湍流模型實現(xiàn)對N-S方程的間接求解，在魚礁型潛堤工程湍流問題中得到廣泛應用的模式是引入雙方程k-ε湍流模型。由于各湍流模型的差異，其結果存在差異，因此根據(jù)試驗結果對引入不同湍流模型的RANS模型和LES模型進行分析比選，以獲得合適的湍流模型和參數(shù)。

2.2.1 基于RANS的雙方程k-ε湍流模型

對于不可壓縮流體，其瞬態(tài)控制方程為

式中：ui為i方向的速度分量；xi為i方向的坐標分量；xj為j方向的坐標分量；t為時間；ρ為流體密度；P為壓強；ν為運動黏滯系數(shù)；fi為作用于單位質(zhì)量水體的體積力。

對瞬態(tài)控制方程進行雷諾平均，得到時均形式的控制方程為

標準k-ε模型是標準的雙方程湍流模型，是由Launder和Spalding［18］通過實驗現(xiàn)象總結出來的半經(jīng)驗公式，其輸運方程可整理如下：

k方程為

ε方程為

式中：k為紊動能；ε為紊動耗散率；Gk為湍流應力或速度梯度引起的k的產(chǎn)生項；Dk為k的耗散項；μ為動力黏滯系數(shù)；μt為渦黏系數(shù)；ρ為密度；σk、σε、Cμ、C1ε、C2ε和C3ε均為經(jīng)驗常數(shù)或改進的變量，具體取值可參考文獻［19］。

標準k-ε模型假定流場為完全湍流，忽略分子之間的黏性，將渦黏系數(shù)假定為各向同性的標量，故其只對完全是湍流的流場有效，并非適用于任何流動。因此，許多學者進一步對其改進，提出不同的假設修正k-ε模型，并由此發(fā)展出許多衍生的k-ε模型。

重整化群k-ε模型是Yakhot和Orszag［20］采用重整化群數(shù)學方法對k-ε模型進行改進后的湍流模型，通過在大尺度運動和修正的黏度項體現(xiàn)小尺度的影響，使這些小尺度運動有系統(tǒng)地從控制方程中去除。重整化群k-ε模型的k方程與標準k-ε模方程中的k方程一樣，但ε方程增加了一項Rε，其源項為計算方法為

同時，在ε方程產(chǎn)生項的系數(shù)C*2ε的計算中引入了主流應變率，這樣在重整化群k-ε模型中C*2ε不僅與流動情況有關，而且在同一問題中也是空間坐標的函數(shù)。

可實現(xiàn)k-ε模型是20世紀末由Shih等［21］提出的帶旋流修正的k-ε模型，與標準k-ε模型相比有2個主要的不同點：一是在模型中為湍流黏性增加了一個公式，湍流黏性系數(shù)仍由式（5）確定，但是Cμ不再是常數(shù)，而是與k、ε和旋轉應變率U*有關的函數(shù)，即

二是耗散率ε輸運方程不同，為耗散率ε增加了新的傳輸項，這個方程來源于一個為層流速度波動而作的精確方程，更適于表達能譜的傳輸［22］?？蓪崿F(xiàn)k-ε模型較前2種k-ε模型的優(yōu)點是可以保持雷諾應力與真實湍流一致，對渦旋流計算、帶壓強梯度的邊界層計算和分離流計算等問題的處理可以更符合實際情況。

2.2.2 大渦模擬

湍流由不同尺度的渦旋組成，湍流的脈動與混合主要由大尺度渦造成，小尺度渦的主要作用是耗散紊動能。大渦模擬以特定的分辨尺度將流動分解為大尺度運動（大于網(wǎng)格尺度）和小尺度運動，對大尺度湍流運動直接用N-S方程求解，而對于小尺度湍流運動通過構建亞格子模型封閉于求解方程組中。不同尺度流動的區(qū)分依靠濾波函數(shù)來實現(xiàn)，對不可壓縮N-S方程在物理空間過濾得到的控制方程為

式中：xi、xj為坐標分量為濾波后的流速分量為濾波后的壓強；τij為亞格子尺度雷諾應力，代表過濾掉的小尺度脈動和大尺度脈動之間的能量運輸。為了使方程閉合，采用Boltzmann假定的渦黏模型求解亞格子應力，為

3 結果分析

3.1 流場分布

圖3為利用PIV技術測量得到的單體魚礁型潛堤近區(qū)在垂向上和水平向上3個不同位置剖面的時均流場圖，通過在流場圖中加入流線的方式更清晰地展示魚礁型潛堤近區(qū)的紊流形態(tài)。由圖3a可以看出，由于礁體的阻水作用，水流行至魚礁型潛堤前方沿著礁體迎流面抬升，在礁體上方形成了明顯的上升流，最大流速達到0.106m·s-1，較入流速度增加了18%。上升流從魚礁型潛堤頂板流過，頂板阻力沿程消耗水流的能量，導致貼近頂板的流速較小。流過魚礁型潛堤頂板的水流，因礁體后沒有了礁體占有的固體空間，水流向下跌落，流線向下，形成了較大的逆時針回旋流，即背渦流，并從背流面底部再次進入礁體。顯然，魚礁型潛堤內(nèi)部水流主要來自從魚礁型潛堤后方進入的背渦流，礁體內(nèi)平均流速較入流速度減小約70%。魚礁型潛堤迎流面和頂板都有透水圓孔，部分水流透過迎流面圓孔進入魚礁型潛堤，產(chǎn)生了類射流效應，在礁體內(nèi)部形成了2個條狀的較高流速區(qū)（vxz≥ 0.058 m·s-1），該區(qū)內(nèi)射流方向與主流方向一致。背渦流底流繞過礁體內(nèi)方柱流向迎流面，受自迎流面圓孔進入的水流干擾，被分割形成上下2對渦旋。同時，由于礁體內(nèi)部流速小，魚礁型潛堤頂部上方流速大，在礁體頂板附近產(chǎn)生了虹吸效應，導致部分水流從頂板圓孔流出，使魚礁型潛堤內(nèi)部的水體與礁頂上方的主流發(fā)生交換。

圖3 水槽試驗結果及特征斷面位置示意Fig.3 Experimental results and the locations of specific cross-section for validation of the numerical models

圖3b和圖3c分別為穿過迎流面近底部（Z=48mm）和近表面（Z=136mm）2排圓孔的水平流場。水流經(jīng)過魚礁型潛堤時，主流過水斷面束窄，在礁體兩側分離形成高流速帶（vxy≥ 0.116 m·s-1）流向下游，并在礁體后方背渦流區(qū)域內(nèi)形成了一對水平展向渦。由于孔隙透水，在魚礁型潛堤內(nèi)部也形成了大小不一的渦旋。在魚礁型潛堤迎流面后方靠近兩側邊界處，受水槽兩側高速水流的影響，礁體內(nèi)、外流體的強力剪切作用最顯著，流速沿魚礁型潛堤兩側邊界的法向變化劇烈，因此在近底部和近表層流場中都出現(xiàn)渦旋，且該處的渦旋結構在不同高度的水平截面上基本一致。近表層流場比近底部流場形成更多渦旋，這是因為近表層截面與魚礁型潛堤頂板的距離小于近底部截面與水槽底部的距離，近表層水流受礁頂壁面影響更大，消耗了更多水流能量。同樣地，在垂向流場中，礁體內(nèi)上半部也產(chǎn)生了更多渦旋，與水平流場情況一致。

為考察不同湍流模型對流場中速度矢量計算結果的精確性，將對流場中不同位置的特征剖面中的縱向流速進行比較。基于PIV試驗結果，提取魚礁型潛堤近區(qū)不同特征斷面沿程的縱向流速，垂向上為穿過礁頂圓孔中心線和背流面的A1—A1、A2—A2和A3—A3；橫向上為穿過礁體前后隔間的中心線C1—C1和C2—C2；縱向上為穿過魚礁型潛堤迎流面圓孔的中心線S1—S1和S2—S2，具體位置如圖3所示。圖4對比了實驗測量值與4種湍流模型的模擬結果。圖4a、4b、4c為3個垂向截面的縱向流速分布圖，其中位于魚礁型潛堤內(nèi)部的A1—A1和A2—A2截面的流速峰值都出現(xiàn)在迎流面透水孔正后方，這也驗證了由于孔隙透水在礁體內(nèi)部產(chǎn)生的射流效應，導致該處水流紊動最為強烈。在A1—A1斷面，各湍流模型對縱向流速變化趨勢的預測與試驗測量結果基本一致，根據(jù)斷面平均流速，與實測值吻合程度的排序依次為標準k-ε＞ 大渦模擬＞ 重整化群kε＞ 可實現(xiàn)k-ε。在A2—A2斷面，標準k-ε和可實現(xiàn)k-ε模型較準確地反映了沿程流速變化的大小和趨勢，可實現(xiàn)k-ε稍優(yōu)于標準k-ε模型，但重整化群k-ε和大渦模擬得到透水孔后方的縱向流速與實測值偏差較大，且遠高于另外2個模型。在A3—A3斷面，由于不再受到迎流面透水孔入流的影響，各湍流模型計算結果均與實測值吻合較好。圖4d、4e為近底層水平流場2個位于中橫向斷面的縱向流速分布圖，各湍流模型對C1—C1斷面流速的計算結果顯示出了較好的吻合性。對于C2—C2斷面，大渦模擬的計算值出現(xiàn)了最大的偏差，大渦模擬得到的C2—C2斷面流速變化與C1—C1的基本一致，在迎流面透水孔正后方都出現(xiàn)了速度峰值，而在其余模型結果中該處流速已大幅降低，無明顯峰值。這由兩方面原因造成，主要原因是大渦模擬水流通過迎流面透水孔后在魚礁型潛堤內(nèi)部形成的條狀高流速區(qū)較其他模型結果長，可延伸至礁體內(nèi)后半?yún)^(qū)，與C2—C2過流斷面相交，導致C2—C2斷面流速在相應位置出現(xiàn)了速度峰值，這在對流場做進一步分析后得到了證實（圖5）；其次，C2—C2斷面位置恰處于2個大小相近、方向相反的渦對交界處（圖3b），由于LES模型對出流斷面流速的變化最敏感［23］，因此在該斷面下大渦模擬計算的流速變化較其他模型大。在其余的3個模型中，以可實現(xiàn)k-ε模型的計算值與實測值最為接近。圖4f、4g為近表層流場中2個穿過迎流面圓孔中心線的縱向流速分布圖，需要注意的是，在魚礁型潛堤迎流面圓孔內(nèi)的流速在試驗測量和數(shù)值計算之間存在較大的差異，產(chǎn)生該現(xiàn)象的原因可能是因為在數(shù)值模型中，壁面小尺度孔結構邊界的準確刻畫對網(wǎng)格模型的要求很高，且近壁面由于邊界層的存在流速變化更為劇烈，即使已針對魚礁型潛堤模型的結構特征做了局部加密處理，但基于計算時間和計算精度綜合考量下選取的網(wǎng)格模型可能仍難以將小尺度孔結構內(nèi)部的實際流速情況精準復演。上述現(xiàn)象產(chǎn)生的原因涉及湍流模型中壁面函數(shù)的選擇、邊界處三維網(wǎng)格劃分以及近壁面處理等多個方面，本文主要討論不同湍流模型對魚礁型潛堤內(nèi)部及近區(qū)復雜湍流情況的預測能力，對這類“結構中的結構內(nèi)部”，即相關魚礁型潛堤中小尺度孔結構內(nèi)流速變化的討論分析仍存在一定局限性，未來需要對其做進一步研究予以解決。除去孔內(nèi)流速外，根據(jù)平均縱向流速分布情況，對于S1—S1，各模型在縱向截面的吻合度由大到小為：標準k-ε、可實現(xiàn)k-ε、 大渦模擬、重整化群k-ε；對于S2—S2，為：可實現(xiàn)k-ε、標準k-ε、重整化群k-ε、大渦模擬。

圖5 魚礁型潛堤近區(qū)流場數(shù)值模擬結果Fig.5 Computed results of flow field in vicinity of the reef-type breakwater

3.2 渦旋結構分布

圖5給出了4個湍流模型對魚礁型潛堤流場的模擬結果，并局部放大了礁體內(nèi)部的流場及其流線變化，以展示各湍流模型對魚礁型潛堤內(nèi)部水流結構的刻畫結果。各模型均較好地預測了魚礁型潛堤前方的上升流和背渦流的形態(tài)，這兩者是最常用于評價魚礁型潛堤流場效應的指標。當垂直方向的流速等于或者大于10%的來流速率時，就把該區(qū)域定義為上升流區(qū)域［6］。表2展示了各湍流模型的計算結果及計算效率，對比垂向上的最大流速可以看出，4種模型的差異很小，并且十分接近實驗數(shù)據(jù)。除大渦模擬模型外，其余模型的相對誤差均小于5%，其中，標準k-ε模型的計算結果誤差最小。流體流過魚礁型潛堤在其背流面形成規(guī)模最大的回流區(qū)，回流區(qū)內(nèi)流向速率為負值，回流區(qū)后流體速率不斷恢復至明渠流分布狀態(tài)。不同湍流模型所得回流區(qū)長度不同，從圖5a可以看出，各模型計算得到的垂向剖面流場圖中魚礁型潛堤后方回流區(qū)的長度各不相同。大渦模擬所得回流長度最短，而重整化群k-ε模型所得回流區(qū)最長，這與陳善群和王澤［24］的研究結論一致。受PIV試驗中CCD相機分辨率的限制，本試驗觀測范圍約為400×300mm，由于魚礁型潛堤后方形成的回流區(qū)規(guī)模較大，本試驗未能捕捉到完整的垂向回流區(qū)，但試驗結果清晰顯示了礁體內(nèi)部的渦流結構特征和礁體后方形成的一對水平展向渦的位置。對比魚礁型潛堤內(nèi)部的垂向渦旋結構，重整化群k-ε和大渦模擬結果較試驗結果偏差較大，主要差別表現(xiàn)在礁體內(nèi)部的前半?yún)^(qū)內(nèi)（迎流面中段正后方），該處流線雖出現(xiàn)較大變形，但仍無渦旋形成，說明上述模型對垂向分離螺旋點的預測較低，這與夏超等［25］的研究結果一致。其余2個湍流模型得到的渦旋結構分布與PIV試驗結果基本吻合，但所得的渦旋數(shù)量存在不同。上述4個湍流模型中，可實現(xiàn)k-ε模型的計算結果在垂向渦旋分布和數(shù)量上都與實測結果最吻合。魚礁型潛堤內(nèi)部的水流結構是由高速入射流和低速緩流相互作用形成的，以連續(xù)的強旋流為主要形態(tài)特征，可實現(xiàn)k-ε模型引入了與旋轉和曲率相關的內(nèi)容，因此在渦旋結構特征明顯的魚礁型潛堤內(nèi)部流場計算中體現(xiàn)了較高的吻合度。

表2 各湍流模型計算結果及計算效率Tab.2 Simulated resulted of different turbulent models and their computational efficiencies

圖5b和圖5c分別為各湍流模型對近底部和近表層水平流場的模擬結果，魚礁型潛堤內(nèi)部和后方渦旋結構的規(guī)模和位置存在一定差異。其中，魚礁型潛堤后方形成的一對展向渦基本以礁體縱向中軸線為對稱。因此以魚礁型潛堤縱向中軸線為基準，計算出不同湍流模型預測的礁體后方右渦核與中軸線和背流面的距離，并通過礁體外部尺寸將其量綱一化，以此來表征渦核的位置（如圖3b所示），所得結果列于表2中。從表中可以看出，在近底層流場中，各湍流模型對渦核位置的預測存在不同程度的誤差，其中，大渦模擬的結果誤差最小，而重整化群k-ε的結果誤差最大，各模型模擬得到的渦核位置距離魚礁型潛堤背流面比實測的更長，這可能是因為試驗中近底部流速較低，部分示蹤粒子沉積在礁體底部，影響了近底層流場結構；而在近表層流場中，標準k-ε模擬結果與試驗結果最接近，大渦模擬次之，重整化群k-ε預測效果仍最差，分析原因為重整化群k-ε中的渦黏系數(shù)μt計算公式中的系數(shù)Cμ是經(jīng)驗常數(shù)，可能只適用于某些流場。進一步對比礁體內(nèi)部的水平渦旋結構，大渦模擬能夠預測最全面的渦旋分布情況，重整化群k-ε亦能模擬得到礁體內(nèi)部較多的大尺度渦旋，但流場內(nèi)的渦旋分布因礁體結構而應呈現(xiàn)的對稱性仍存在偏差；可實現(xiàn)k-ε模型對近底部水平流場的預測效果優(yōu)于近表層，但在近表層水平流場中，僅該模型復演了與PIV測試結果一致的附著于礁體背流面方柱的2個尺度最小的角渦；標準k-ε的模擬結果誤差最大，在礁體內(nèi)部預測到的渦旋數(shù)量最少，在近表層中甚至表現(xiàn)為無渦旋流場。這是由于標準k-ε模型假定了渦黏系數(shù)μt是各向同性的標量，所以其用于模擬魚礁型潛堤內(nèi)部復雜結構引起的強渦旋流動時會產(chǎn)生一定的失真。

由表2 中可見，模型計算效率大小次序為：可實現(xiàn)k-ε、標準k-ε、重整化群k-ε、大渦模擬，由于大渦模擬模型對網(wǎng)格精度和計算資源的要求最高，故計算效率最低。雖然大渦模擬計算量相對于k-ε模型大，但大渦模擬方法對N-S方程直接求解大渦的優(yōu)勢使其對礁體近區(qū)渦旋結構的預測效果明顯優(yōu)于k-ε模型，尤其是魚礁型潛堤內(nèi)部和后方的展向渦對形態(tài)。雙方程k-ε模型對上升流的模擬效果比背渦流更加準確，其中雙方程模型以可實現(xiàn)k-ε模型對礁體近區(qū)的流動分離現(xiàn)象預測最準確，且計算效率最高?？傮w而言，在計算能力有限的條件下，可優(yōu)先考慮使用可實現(xiàn)k-ε模型對多孔箱體魚礁型潛堤的三維湍流場進行模擬計算。

4 結論

針對多孔箱體魚礁型潛堤，分別采用水槽試驗和數(shù)值模擬獲得了單體魚礁型潛堤近區(qū)的三維繞流流場和渦旋結構，結果表明礁體近區(qū)復雜的繞流流場是由不同強度的渦旋之間相互作用的結果。此外，基于試驗結果，評估了4種典型的湍流模型的計算精度和效率，所得主要結論如下：

（1）試驗結果表明，魚礁型潛堤前形成了明顯的上升流，最大上升流速為0.106 m·s-1，較入流速率增加了18%；背渦流在礁體后方形成，水平上呈現(xiàn)為一對以魚礁型潛堤縱向中軸線為對稱的展向渦對，并穿過背流面進入魚礁型潛堤，礁內(nèi)平均流速較入流速率減小約70%；水流通過迎流面圓孔以射流方式進入魚礁型潛堤，在對應圓孔后方形成了2個較高流速區(qū)（vxz≥ 0.058 m·s-1），并對礁內(nèi)水流結構起到重新調(diào)整的作用，在垂向上將內(nèi)部水體分割成了上下2對渦旋，近表層由于受礁頂壁面影響更大而產(chǎn)生了更多渦旋。

（2）采用的4種湍流模型都能在一定精度內(nèi)模擬魚礁型潛堤繞流流場的速度分布，其中標準k-ε和可實現(xiàn)k-ε在縱向流速的計算上精度最高；對魚礁型潛堤上升流的預測各模型差別不大，且都接近實測值；對于背渦流的復演，各湍流模型對渦核位置的預測都存在不同程度的誤差，大渦模擬誤差最小，而重整化群k-ε所得流場偏差最大；此外，雙方程k-ε模型對上升流的模擬效果比背渦流更準確。

（3）礁體內(nèi)部形成的渦旋結構，各湍流模型之間呈現(xiàn)出較大的差別。對于垂向渦旋結構，重整化群k-ε和大渦模擬對垂向分離螺旋點的預測較低，可實現(xiàn)k-ε模型的計算結果在垂向渦旋分布和數(shù)量上都與實測結果最吻合，且模型計算效率最高；對于水平渦旋結構，大渦模擬能夠預測最全面的渦旋分布，但模型計算效率最低，標準k-ε在紊流場的刻畫效果最低，可實現(xiàn)k-ε對近壁面小尺度角渦的刻畫效果最佳。

（4）由于透水魚礁型潛堤結構的復雜性，考慮到礁體近區(qū)繞流流動是三維的，且流態(tài)復雜，礁體內(nèi)部流速很低，準確模擬存在難度。在比選的4個湍流模型中，綜合模型的計算精度和效率來看，可實現(xiàn)k-ε是在計算能力有限的條件下最適用于多孔箱體魚礁型潛堤近區(qū)復雜繞流流場模擬的湍流模型。

作者貢獻聲明：

匡翠萍：研究思路與論文結構確定、論文審閱和修改。

鄭宇華：實驗操作、數(shù)據(jù)處理與分析、初稿撰寫與修改。

顧 杰：實驗設計與指導。

韓雪?。簩嶒灢僮鳌?/p>