王 寅，王玲玲，計(jì) 勇，席芝橙，張?chǎng)?/p>

(1.南昌工程學(xué)院水利與生態(tài)工程學(xué)院，江西 南昌 330099; 2.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

內(nèi)波通常發(fā)生于密度穩(wěn)定分層的水體內(nèi)部[1]。在海洋、河口、湖泊等流體系統(tǒng)中，溫度、鹽度以及其他環(huán)境因素引起的密度沿水深方向上的改變，會(huì)導(dǎo)致水體密度的分層[2]。在此類(lèi)環(huán)境中，微弱的外界擾動(dòng)就可能在密度分界面(密度躍層)上引發(fā)攜帶巨大能量的內(nèi)波現(xiàn)象。內(nèi)孤立波是一種特殊的非線(xiàn)性?xún)?nèi)波，其振幅大、周期短、所攜能量大[3]，大振幅內(nèi)孤立波對(duì)近岸及河口復(fù)雜水動(dòng)力環(huán)境中工程結(jié)構(gòu)物的安全穩(wěn)定性有非常大的影響[4]。內(nèi)孤立波在傳播過(guò)程中以密度躍層為界上下水層呈反向流動(dòng)，由此產(chǎn)生的剪切流所攜帶的巨大能量會(huì)導(dǎo)致異常強(qiáng)大的流速。內(nèi)孤立波2.1m/s波致流作用下的柱體所承受的最大荷載相當(dāng)于波長(zhǎng)300 m、波高達(dá)到18 m的表面波作用[5]，這種強(qiáng)烈的作用力可能對(duì)近岸及河口復(fù)雜水動(dòng)力環(huán)境中水下立管及水下支撐樁柱的安全穩(wěn)定性造成嚴(yán)重的威脅[6]。深入研究?jī)?nèi)孤立波對(duì)水下結(jié)構(gòu)物的作用機(jī)理以及結(jié)構(gòu)物內(nèi)波動(dòng)力學(xué)問(wèn)題具有非常重要的工程應(yīng)用價(jià)值。

柱體所受水平力可分為壓差阻力和摩擦阻力兩部分[7]。壓差阻力是由于柱體迎流面和背流面的壓差所產(chǎn)生的水平力；摩擦阻力是由于流體黏性作用所產(chǎn)生的水平力，即黏滯力。以往受力研究多只計(jì)算出內(nèi)孤立波對(duì)柱體的整體作用力[8-10]，并未考慮在不同分層水體中柱體不同部分受力隨水深分布情況，因此并不能獲知柱體詳盡的受力特性。本文借助三維數(shù)值波浪水槽，采用大渦模擬(large-eddy simulation，LES)技術(shù)模擬內(nèi)孤立波的產(chǎn)生和傳播，通過(guò)提取數(shù)模計(jì)算結(jié)果中不同水層柱體表層的壓強(qiáng)分布及切向速度，采用微元法積分出柱體在各水層所受壓差阻力以及黏滯力，由此獲取受力沿水深分布特性，以揭示分層流環(huán)境下柱體的受力機(jī)理。

1 控制方程及求解方法

1.1 動(dòng)量輸運(yùn)方程

基于連續(xù)性假設(shè)，可用Navier-Stokes方程(N-S方程)描述不可壓縮黏性流體三維瞬態(tài)運(yùn)動(dòng)過(guò)程：

(1)

(i，j=1，2，3)

(2)

式中：ρ為密度項(xiàng)；t為時(shí)間；xi、xj為笛卡爾坐標(biāo)系的3個(gè)方向坐標(biāo)；ui、uj為笛卡爾坐標(biāo)系中的3個(gè)流速分量；p為壓力項(xiàng)；μ為動(dòng)力黏性系數(shù)；fi為i方向的單位體積力。

1.2 標(biāo)量輸運(yùn)方程

采用考慮了兩層水體之間的質(zhì)量交換作用以及對(duì)流-擴(kuò)散影響的標(biāo)量輸運(yùn)方程：

(3)

ρ=C+ρ0

(4)

式中：C為鹽度，kg/m3；k為分子擴(kuò)散系數(shù)；S為源項(xiàng)；ρ0為清水密度，kg/m3。

1.3 大渦模擬

大渦模擬技術(shù)是對(duì)紊流脈動(dòng)的一種空間平均，通過(guò)濾波函數(shù)將大尺度的渦和小尺度的渦分離開(kāi)[11-14]。過(guò)濾后的動(dòng)量和質(zhì)量方程可表示為

(5)

(6)

(7)

1.4 柱體壓差阻力的求解方法

1.4.1圓柱

如圖1所示，P為水深為y的截面圓柱表面上受到的點(diǎn)壓強(qiáng)(基于數(shù)值模擬結(jié)果提取)，則該截面圓周所受水平方向的壓差阻力fc可定義為

(8)

式中：θ為點(diǎn)壓強(qiáng)P與x軸的夾角，(°)，以逆時(shí)針為正；l為弧長(zhǎng)，cm。

圖1 圓柱壓差阻力計(jì)算方法示意圖

基于微元法，可將式(8)離散為

(9)

式中：Pi為作用在第i個(gè)微元上的壓力，Pa；θi為第i個(gè)微元法線(xiàn)延長(zhǎng)線(xiàn)與x軸的夾角，(°)，以逆時(shí)針為正；n為圓周被等分?jǐn)?shù)，本文取n=360；Δli為第i個(gè)微元的弧長(zhǎng),Δli=2πr/n=0.087 cm(r為圓柱半徑，本文取r=5 cm)。

1.4.2方柱

如圖2所示，假定水深為y的截面方柱表面上受到的點(diǎn)壓強(qiáng)為P，則該方柱截面所受水平壓差阻力fs定義為

(10)

式中：fsj為作用在方柱截面第j(j=1，2，3，4)條邊上的壓力，N；Pi為方柱各邊第i個(gè)微元上的壓力，N；Δsi為第i個(gè)微元的邊長(zhǎng)，Δsi=4s/m=0.11 cm(s為方柱周長(zhǎng)，本文取s=10 cm;m為方柱周長(zhǎng)被等分?jǐn)?shù)，本文取m=360)。

圖2 方柱壓差阻力的計(jì)算方法示意圖

1.5 柱體黏滯力的求解方法

1.5.1圓柱

如圖3所示，可將距離圓柱表面最近一層網(wǎng)格上某微元a的切向速度uτ和徑向速度uγ(基于數(shù)值模擬結(jié)果提取)定義為

uτ=uzcosθ-uxsinθ

(11)

uγ=uzsinθ+uxcosθ

(12)

式中：ux、uz分別為該微元的x方向和z方向速度，cm/s；θ為微元a與x軸的夾角，(°)，以逆時(shí)針為正。則某微元i所受水平方向上的黏滯力Δfcτi為

(13)

圖3 圓柱黏滯力計(jì)算方法示意圖

基于微元法，圓柱表面所受水平方向上的黏滯力fcτ為

(14)

1.5.2方柱

如圖4所示，沿用上述圓柱黏滯力的計(jì)算方法，可求得水深為y的截面方柱表面所受水平方向上的黏滯力fsτ：

圖4 方柱黏滯力計(jì)算方法示意圖

(15)

其中

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

建立三維數(shù)值波浪水槽如圖5所示，水槽長(zhǎng)12 m、寬0.6 m、高0.8 m，柱體周?chē)骄W(wǎng)格和水槽垂向網(wǎng)格劃分方法如圖6、圖7所示。

圖5 三維數(shù)值波浪水槽示意圖(單位：m)

圖6 柱體周?chē)骄W(wǎng)格劃分方法

圖7 垂向網(wǎng)格劃分方法

定義在數(shù)值模擬中柱體所受的無(wú)量綱水平合力F′的計(jì)算公式[15]為

(16)

圖8 不同波幅時(shí)柱體的F′隨時(shí)間t的變化過(guò)程

式中：F為數(shù)值模擬中計(jì)算所得的柱體所受水平力合力，N；A為柱體的迎風(fēng)面積，cm2；H為水槽總水深，cm；g為重力加速度，m/s2。

2.1 計(jì)算工況

數(shù)值模擬采用的波形為下凹型，參數(shù)設(shè)定及工況選定具體為：上層為清水，密度ρ1=0.998 g/cm3，水層厚度h1=0.2 m；下層為鹽水，密度ρ2=1.017 g/cm3，水層厚度h2=0.6 m，總水深H=0.8 m，設(shè)定4種工況N1、N2、N3、N4，相對(duì)波幅ηo/H分別為0.026、0.044、0.054和0.098。如圖9所示，采用重力塌陷法制造內(nèi)孤立波[16]，圖中造波區(qū)兩側(cè)密度躍層高度差ho稱(chēng)為階躍高度，造波區(qū)長(zhǎng)度Lo稱(chēng)為階躍長(zhǎng)度。

圖9 重力塌陷法造波示意圖

2.2 數(shù)值模擬方法和邊界條件

三維數(shù)值波浪水槽中下凹型內(nèi)孤立波的形成和傳播采用大渦模擬技術(shù)模擬；利用SIMPLE算法對(duì)流速-壓力項(xiàng)進(jìn)行耦合，以保證質(zhì)量守恒，并由此得到壓力場(chǎng)；利用二階中心差分格式及隱式格式分別對(duì)空間和時(shí)間進(jìn)行離散。造波區(qū)端邊(圖9中的左邊壁)、數(shù)值水槽底部以及柱體壁面均采用無(wú)滑移固壁邊界，水槽兩側(cè)采用滑移固壁邊界。為了使內(nèi)波不發(fā)生反射，水槽右端采用Sommerfeld輻射型邊界[17]；頂部采用“剛蓋”假定，以此忽略表面波的影響[17]。

2.3 數(shù)學(xué)模型驗(yàn)證

采用與文獻(xiàn)[18-19]相同的方法驗(yàn)證了建立的數(shù)學(xué)模型，比較了不同波幅情況下，圓柱和方柱模型試驗(yàn)與數(shù)值模擬中內(nèi)孤立波作用力的時(shí)歷曲線(xiàn)。結(jié)果表明模型試驗(yàn)與數(shù)值模擬結(jié)果吻合較好，驗(yàn)證了建立的數(shù)學(xué)模型的可靠性，可用于后續(xù)試驗(yàn)的模擬。

3 模擬結(jié)果與分析

內(nèi)孤立波環(huán)境下的速度場(chǎng)和壓力場(chǎng)與密度均一流環(huán)境區(qū)別較大，在內(nèi)孤立波作用下柱體受力三維特性更加明顯[18]。準(zhǔn)確剖析柱周流場(chǎng)及柱表壓強(qiáng)分布是科學(xué)研究和預(yù)測(cè)近岸及河口地區(qū)樁柱受力的必要條件。柱體所受水平力可分為壓差阻力和摩擦阻力(黏滯力)兩部分，這兩類(lèi)力求解的基本思路：①在柱身沿垂直(y)方向每隔一定距離截取1個(gè)斷面，基于模擬結(jié)果在每個(gè)斷面上沿柱周向選取360個(gè)點(diǎn)，并提取這些點(diǎn)的壓強(qiáng)值及流速值(得到柱周的壓強(qiáng)分布及流速分布)；②采用微元法，通過(guò)積分不同斷面柱周的壓強(qiáng)分布和流速分布可分別得到圓柱及方柱的壓差阻力f(式(9)(10))和黏滯力fτ(式(14)(15))。

3.1 柱體壓差阻力沿水深的分布

定義柱體在水深y處所受無(wú)量綱水平壓差阻力f′為

(17)

3.2 柱體黏滯力沿水深的分布

(18)

表1 圓柱和方柱的無(wú)量綱化水平壓差阻力

3.3 單柱受力的量化分析

表2 柱體受力量化分析工況設(shè)置及相關(guān)參數(shù)

(19)

m=-0.88h1/h2+3.72η0/H+0.34

(20)

(21)

此時(shí)R2=0.957，擬合優(yōu)度高(圖10)。綜合影響系數(shù)m雖是擬合過(guò)程中得到的，但其具有一定的物理意義：內(nèi)孤立波相對(duì)波幅越大、上下水層厚度差距越大，則m值越大，對(duì)應(yīng)柱體受力越大。m為兩參數(shù)η0/H和h1/h2對(duì)柱體受力的影響提供了量化的指標(biāo)，且m與η0/H正相關(guān)，與h1/h2負(fù)相關(guān)(式(20))。綜上可知，柱體受力在不同分層條件下均隨相對(duì)波幅η0/H的增大而增大；而在同一波幅條件下，水深比h1/h2越大，柱體受力越小。

圖與ηo/H和h1/h2的相關(guān)關(guān)系

4 結(jié) 論

a. 波幅相同時(shí)，以密度躍層為界，方柱于上下水層中所受到的壓差阻力均大于圓柱：上層水體中方柱受力約為圓柱的1.5倍；下層水體中方柱受力約為圓柱的3.5倍。

b. 作用于柱體上的黏滯力比壓差阻力小1或2個(gè)數(shù)量級(jí)，因此可忽略流體的黏性效應(yīng)，只考慮壓差阻力的作用效果。

c. 受力峰值系數(shù)與相對(duì)波幅及水深比具有較強(qiáng)的相關(guān)性，且與綜合影響參數(shù)呈明顯的指數(shù)分布關(guān)系，兩者擬合優(yōu)度為0.957，擬合優(yōu)度高。受力峰值系數(shù)與相對(duì)波幅正相關(guān)，與水深比負(fù)相關(guān)。