二維造流水槽整流裝置數(shù)值模擬

馬 暄，白志剛

（天津大學(xué) 建筑工程學(xué)院，天津 300350）

1 研究背景

波流水槽是港口、水利和海洋船舶工程及相關(guān)學(xué)科的重要實(shí)驗(yàn)設(shè)備。為兼具造流、造波的功能，常見的波流水槽一般采用上部造波，下部造流的結(jié)構(gòu)形式。由于造流廊道出入口處流態(tài)復(fù)雜，可供用于物模實(shí)驗(yàn)的理想?yún)^(qū)域僅限于水槽中部，有效試驗(yàn)段長(zhǎng)度只有50%左右。為了盡可能增加試驗(yàn)段，實(shí)驗(yàn)室建設(shè)中常采用加大水槽長(zhǎng)度和縮短實(shí)驗(yàn)區(qū)長(zhǎng)度的方法保證實(shí)驗(yàn)結(jié)果的有效性，但前者加大了建設(shè)投資，后者則降低了實(shí)驗(yàn)精度。為了改善水槽出入口處的流態(tài)，實(shí)驗(yàn)中常采用盲溝材料作為整流裝置。Tollmien等［1］于1961年首次提出了采用整流網(wǎng)改善流場(chǎng)品質(zhì)的方法。Corrsin［2］在物理模型實(shí)驗(yàn)中驗(yàn)證了Tollmien的設(shè)想。Laws E M［3］等以金屬絲-紗網(wǎng)為例，研究了整流網(wǎng)對(duì)時(shí)均流速和湍流分布的影響。

盲溝材料雖然廣泛應(yīng)用于實(shí)驗(yàn)水槽，但盲溝材料的孔隙率，整流段的長(zhǎng)度等參數(shù)并沒有明確的標(biāo)準(zhǔn)，實(shí)際建設(shè)中常依賴經(jīng)驗(yàn)確定。建立相應(yīng)的數(shù)學(xué)模型有助于確定盲溝整流裝置對(duì)流場(chǎng)的作用效果。由于盲溝材料結(jié)構(gòu)復(fù)雜，建模難度大。且水流流經(jīng)盲溝材料時(shí)流態(tài)較為復(fù)雜，流態(tài)分布難以模擬。對(duì)于盲溝材料的數(shù)學(xué)模型的研究尚不充分，具有一定的研究?jī)r(jià)值。

SPH方法在大變形流體問題的模擬當(dāng)中具有很好的自適應(yīng)性［4］。在流速邊界模擬方面，本文采用由Federico等［5］提出的周期性邊界法，通過將水槽末端的粒子轉(zhuǎn)移到水槽前端，實(shí)現(xiàn)粒子的循環(huán)流動(dòng)。該方法可以保證模擬區(qū)域內(nèi)粒子總數(shù)恒定不變，模擬過程中不需要?jiǎng)?chuàng)建新的水粒子。在阻尼區(qū)域的處理方面，數(shù)值阻尼法是SPH方法中常用的邊界處理方法。其主要原理是通過在部分區(qū)域逐漸增加黏性，從而降低質(zhì)點(diǎn)速度。該方法原本主要用于消除波浪模擬中的反射波［6］。本文通過修改數(shù)值阻尼區(qū)域，將其應(yīng)用在造流模擬中?；贠wen等的整流理論，用數(shù)值阻尼區(qū)域的宏觀特性表征實(shí)際篩網(wǎng)材料的微觀特性，調(diào)節(jié)水流流速分布，達(dá)到實(shí)際整流裝置的效果。通過與物理模型實(shí)驗(yàn)相對(duì)比，數(shù)值阻尼區(qū)域的整流效果與物理模型實(shí)驗(yàn)結(jié)果偏差在10%以內(nèi)，具有較好的真實(shí)性和有效性。

2 水槽數(shù)值模擬方法

2.1 SPH方法簡(jiǎn)介 SPH是一種無(wú)網(wǎng)格的拉格朗日方法，在大變形流體和自由表面流體的數(shù)值模擬方面優(yōu)于歐拉方法。SPH方法的基本方程是拉格朗日形式的動(dòng)量方程：

式中：v為質(zhì)點(diǎn)速度；ρ為質(zhì)點(diǎn)密度；P為作用在質(zhì)點(diǎn)上的壓力的合力；g為重力加速度；Γ為包括黏性力、切應(yīng)力等在內(nèi)的耗散項(xiàng)的總和。針對(duì)不同的問題，可以對(duì)Γ做不同的處理。每個(gè)計(jì)算時(shí)間步中，通過式（1）得到每個(gè)粒子的加速度，從而計(jì)算出粒子的位移。計(jì)算域中所有的粒子都會(huì)參與計(jì)算，對(duì)于流體粒子，計(jì)算出的位移即為下一個(gè)時(shí)間步粒子的新位置；而對(duì)于邊界粒子，計(jì)算出的位移會(huì)在時(shí)間步內(nèi)清零，即視為剛性邊界。

粒子密度ρ通過式（2）計(jì)算。對(duì)于粒子a來(lái)說(shuō)，近似認(rèn)為該粒子只受到以該粒子為圓心，半徑h以內(nèi)的所有其他粒子的影響，h稱為光滑長(zhǎng)度。

式中：t為時(shí)間；mb為質(zhì)點(diǎn)密度；vab為va-vb即粒子速度的矢量差，a和b為粒子編號(hào)；Wab為核函數(shù)。核函數(shù)的選取決定了算法的精度及計(jì)算速度，本文選用了Wendland五次核函數(shù)［7］，其數(shù)學(xué)式為：

用Monaghan J J的方法［8］將式（1）中的壓力項(xiàng)離散化。

SPH方法常用人工黏性表征水體的黏性，但本文中流場(chǎng)雷諾數(shù)較大，考慮到紊動(dòng)切應(yīng)力的存在，將耗散項(xiàng)Γ表示為層流黏性項(xiàng)和紊流黏性項(xiàng)之和。

層流黏性［11］表示為

紊流黏性［12］表示為

將式（4）—式（7）代入式（1）得：

本數(shù)模的模擬時(shí)間長(zhǎng)達(dá)120 s，故采用二階顯式Symplectic算法［13］，可在長(zhǎng)時(shí)間模擬中保持較高的計(jì)算精度。

2.2 水動(dòng)力邊界條件 水槽的固壁邊界采用動(dòng)力學(xué)邊界粒子法［14］，在固壁邊界上排布2層與流體粒子相同屬性的邊界粒子，通過與流體粒子間的連續(xù)性方程求得其物理量。具有計(jì)算上的便捷性和對(duì)復(fù)雜邊界問題的適應(yīng)性。

為使計(jì)算區(qū)域中的流體質(zhì)點(diǎn)保持穩(wěn)定的循環(huán)流動(dòng)，將周期性邊界條件和流速邊界條件相結(jié)合，模擬了物理模型實(shí)驗(yàn)中水泵的功能。首先將流出計(jì)算區(qū)域的粒子從計(jì)算域的另一端重新加入流場(chǎng)，使模擬過程中計(jì)算區(qū)域中的粒子總數(shù)不變。同時(shí)，在周期性邊界邊緣的粒子也會(huì)與另一側(cè)周期性邊界邊緣的粒子相互作用，避免了邊界粒子缺失的問題［15］。在粒子的出入口一定長(zhǎng)度的區(qū)域內(nèi)的粒子強(qiáng)制賦予流速，使得計(jì)算區(qū)域中的粒子受到流速邊界處粒子的影響，而流速邊界處粒子則不受到計(jì)算區(qū)域粒子的影響［4］，如圖1所示。

圖1 邊界粒子分布示意圖

2.3 數(shù)值整流區(qū)域 疏密不均的整流網(wǎng)可將均勻流非均勻化，亦可將非均勻流均勻化。若假設(shè)流體為無(wú)旋流，則整流網(wǎng)的流速重分布作用主要與整流網(wǎng)阻流系數(shù)K和上游流速U有關(guān)。以剪切流為例，若整流網(wǎng)上游為流速為U的均勻流，下游產(chǎn)生流速分布為u的剪切流，如圖2所示。流速u(y)分布規(guī)律為：

圖2 網(wǎng)格分布與流速分布的關(guān)系［16］

系數(shù)λ滿足：

整流網(wǎng)對(duì)流體湍流度具有調(diào)節(jié)作用。整流網(wǎng)既可提高流體的湍流度，亦可降低流體的湍流度。整流網(wǎng)下游的湍流度主要由上游流體透過整流網(wǎng)的湍流、整流網(wǎng)本身產(chǎn)生的湍流和下游流體不勻均流衰減過程中產(chǎn)生的湍流組成［3］。其中整流網(wǎng)產(chǎn)生的湍流主要與整流網(wǎng)的孔隙率ξ和基于絲徑的雷諾數(shù)Rd有關(guān)。對(duì)于盲溝整流網(wǎng)，當(dāng)ξ≤ 0.3且Rd＜40時(shí)，整流網(wǎng)本身不會(huì)產(chǎn)生湍流［17］。

整流裝置前后的速度平方差為：

在計(jì)算域中設(shè)定阻流區(qū)域，當(dāng)粒子進(jìn)入該區(qū)域時(shí)，除粒子間相互作用以外，還會(huì)受到阻流區(qū)域的作用。

式中：系數(shù)k為阻流系數(shù)，l為整流段的厚度；α為修正系數(shù)，通過試算確定；θ為設(shè)定流向與水平方向的夾角，通過設(shè)定合適的θ，可令阻流區(qū)域只約束粒子在指定方向的速度。由于阻流區(qū)域在垂直方向阻流系數(shù)不相同，粒子經(jīng)過阻流區(qū)后即會(huì)產(chǎn)生不同的流速梯度。通過調(diào)節(jié)合適的修正系數(shù)，即可模擬實(shí)際整流網(wǎng)的整流效果。

3 水槽物理模型實(shí)驗(yàn)

為了驗(yàn)證盲溝材料在造流水池中的整流效果，在天津大學(xué)“水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室”建設(shè)如圖3所示的造流水槽物理模型。模型水槽基本尺寸與數(shù)學(xué)模型中的相同，水流由圖3中左側(cè)水泵推動(dòng)之后，通過進(jìn)水口進(jìn)入水槽，再由出水口經(jīng)右側(cè)水泵，通過水槽下方的回水廊道回到左側(cè)?；厮鹊纼?nèi)設(shè)有回流水泵，用于補(bǔ)償水流流經(jīng)回水廊道的能量損失，使a處和f處的水位保持一致。由于水槽的兩端各設(shè)有一個(gè)相同的水泵，可認(rèn)為水槽兩端具有相同的流速邊界條件，與數(shù)學(xué)模型中的流速邊界吻合。

圖3 造流水槽物理模型整體布置圖

實(shí)驗(yàn)水槽中的水泵流量與水泵電機(jī)頻率和水泵兩側(cè)的水位差有關(guān)。通過預(yù)實(shí)驗(yàn)試測(cè)得到水泵流量與水槽內(nèi)水深及電機(jī)頻率的對(duì)應(yīng)關(guān)系如表1和圖4所示。

表1 水泵流量與水深和頻率對(duì)應(yīng)關(guān)系 （單位：m3/s）

圖4 水泵流量對(duì)應(yīng)關(guān)系圖

實(shí)際流入水槽的流量并不等于水泵的造流流量，而是滿足式（13），如圖5所示。

圖5 泵室水流流態(tài)示意圖（單位：Hz）

隨著水泵頻率增大，水泵槳葉旋轉(zhuǎn)變快，廊道抽水流量Q1隨之增大，出流口實(shí)際流量Q2增大。隨著水槽內(nèi)水深增大，水泵內(nèi)部流量略有增大，這是由于水泵的安裝方式屬于濕式安裝，水泵出流口與水槽進(jìn)口之間存在縫隙。水深較小時(shí)，泵口回流流量Q′2增大，出流口實(shí)際流量Q2減小，即實(shí)際流入水槽內(nèi)的流量較小。

將盲溝整流網(wǎng)布設(shè)在水槽的末端，如圖6所示。分別將水深調(diào)整為25、35、45和55 cm，在距離進(jìn)水口2.25、4.25、6.25、8.25和10.25 m處各布置一個(gè)流速傳感器，測(cè)量水下距離水面10 cm處流速，如圖7所示。

圖6 盲溝整流網(wǎng)安裝位置

圖7 流速測(cè)點(diǎn)布置圖（單位：mm）

4 水槽數(shù)值計(jì)算分析

4.1 數(shù)模整體布置及有效性分析 如圖8所示，模型水槽總長(zhǎng)22.19 m，中部平臺(tái)底長(zhǎng)14.62 m，水池深1.25 m。模型有效性主要與水槽入口處的流速有關(guān)。理論上，深度方向粒子布置得越多，計(jì)算結(jié)果越理想。取粒子間距的最不利工況（水深25 cm），先后將粒子間距設(shè)置為0.04、0.03、0.02、0.01和0.005 m分別試算。分別測(cè)量各個(gè)粒子間距條件下圖8中陰影區(qū)域的垂向流速分布，繪制在半對(duì)數(shù)坐標(biāo)系中，如圖9所示。

圖8 造流水槽計(jì)算模型整體布置圖 （單位：m）

如圖9所示，當(dāng)粒子間距小于0.01 m時(shí)，流速分布不再受到粒子間距的影響，可以認(rèn)為0.01是足夠有效的粒子間距?？紤]計(jì)算耗時(shí)，取0.01 m為模型的粒子間距。

圖9 不同粒子間距入口垂向流速分布對(duì)比圖

數(shù)模中的流速分布需與物模實(shí)驗(yàn)結(jié)果相對(duì)應(yīng)。取水深的最不利工況（水深75 cm），將數(shù)模和物模中的入口流速垂向分布繪制在半對(duì)數(shù)坐標(biāo)系中，結(jié)果如圖10所示。由圖10可知，數(shù)模中的入口流速垂向分布與物理模型實(shí)驗(yàn)結(jié)果吻合較好，且整體上呈對(duì)數(shù)分布，符合紊流的流速分布特征?？梢哉J(rèn)為本數(shù)模中的流速分布真實(shí)有效。

圖10 出流口處垂向流速分布對(duì)比圖

4.2 實(shí)驗(yàn)水槽流速分布對(duì)比 測(cè)量有整流裝置時(shí)的水槽內(nèi)部沿程流速分布情況，分別與有無(wú)整流裝置的數(shù)模結(jié)果比對(duì)，如圖11所示。

圖11中，當(dāng)水深小于35 cm時(shí)，無(wú)論是否加入數(shù)值整流區(qū)，數(shù)模中物模的流速數(shù)據(jù)均擬合較好；當(dāng)水深大于45 cm時(shí)，數(shù)值整流區(qū)對(duì)流速分布影響顯著。無(wú)整流區(qū)時(shí)，55 cm水深的數(shù)模物模差異達(dá)到40.36%，有整流區(qū)時(shí)，數(shù)模物模差異在小于10%。這是因?yàn)楦咚顥l件下，流速測(cè)量點(diǎn)剛好位于流體上側(cè)的低流速區(qū)域。如圖12和圖13所示，圖中不同灰度表示不同流速分布，流速隨顏色深度加大?！啊痢北硎揪嚯x進(jìn)水口10250處的測(cè)點(diǎn)位置。無(wú)整流區(qū)的組次中，由于測(cè)點(diǎn)位于上側(cè)的低流速區(qū)流速偏小，表現(xiàn)為圖11（d）中10 250處的最小值；有整流區(qū)的組次中，由于低流速區(qū)長(zhǎng)度得到有效控制，測(cè)量流速與物模結(jié)果吻合。水深55 cm，頻率25 Hz條件下各水平位置流速分布如圖14所示。

圖11 不同水深條件下不同水平位置流速對(duì)比圖

圖13 55cm水深條件下有數(shù)值整流區(qū)域時(shí)水槽末端流速分布著色圖

圖14 55cm水深條件下不同水平位置流速分布圖

5 結(jié)論

（1）將流速邊界條件和周期性邊界條件相結(jié)合，形成了穩(wěn)定的流場(chǎng)。通過與物模實(shí)驗(yàn)中進(jìn)水口的流速分布相對(duì)比，說(shuō)明了數(shù)學(xué)模型的合理性。結(jié)果表明，當(dāng)粒子間距足夠小時(shí)，可以得到穩(wěn)定、合理的流速場(chǎng)。

（2）整流網(wǎng)相關(guān)理論可應(yīng)用在SPH方法當(dāng)中，模擬實(shí)際整流網(wǎng)對(duì)剪切流的整流效果。本文根據(jù)未加整流區(qū)域時(shí)的流速分布規(guī)律設(shè)置了數(shù)值整流區(qū)域，并驗(yàn)證了整流效果。物模實(shí)驗(yàn)顯示，有整流網(wǎng)的工況中水槽中部流速分布均勻。數(shù)模實(shí)驗(yàn)顯示，未加數(shù)值整流區(qū)域時(shí)，流速分布與物模結(jié)果相差較大；加入整流區(qū)域后，流速分布與物模結(jié)果對(duì)應(yīng)較好。

（3）整流網(wǎng)的整流效果與水深有關(guān)。當(dāng)水深小于35 cm時(shí)，不同水深的流速差異較小，數(shù)值整流區(qū)域?qū)α魉俚挠绊懖⒉幻黠@。當(dāng)水深大于45 cm時(shí)，未添加數(shù)值整流區(qū)域的數(shù)值模擬結(jié)果與物理模型實(shí)驗(yàn)差異較大，添加數(shù)值整流區(qū)域的數(shù)值模擬結(jié)果與物理模型實(shí)驗(yàn)結(jié)果相符。因此，本文的數(shù)值整流方法可顯著改善大水深條件下的垂向流速均勻性。

需指出的是，受實(shí)驗(yàn)條件限制，本文未能測(cè)量出物模實(shí)驗(yàn)中無(wú)盲溝整流材料時(shí)水槽尾端的流速分布，而是通過有效性驗(yàn)證的數(shù)模結(jié)果作為替代。且本文未能測(cè)量出物模實(shí)驗(yàn)中水槽末端的流速的局部分布，而是通過水槽內(nèi)部的流態(tài)間接反映水槽尾部的流態(tài)。該缺陷可能導(dǎo)致數(shù)值阻尼區(qū)域的阻流系數(shù)取值存在偏差。該偏差可通過調(diào)整修正系數(shù)α作為彌補(bǔ)。