劉 宣，白志剛，余海濤，高樹(shù)軍

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

波浪水槽是研究波浪傳播及其與結(jié)構(gòu)物直接相互作用的重要手段之一，具有直觀、可靠等優(yōu)點(diǎn)。其缺點(diǎn)是成本高、過(guò)程復(fù)雜同時(shí)還可能存在比尺問(wèn)題。隨著計(jì)算機(jī)技術(shù)的發(fā)展和計(jì)算流體力學(xué)理論的進(jìn)步，基于N-S方程建立的粘性數(shù)值水槽得到了迅速發(fā)展，在模擬波浪破碎、旋渦，波浪與建筑物相互作用等方面問(wèn)題有著很好的應(yīng)用。Miyata[1]等采用MAC法追蹤自由面，有限差分法離散N-S方程，建立二維數(shù)值波浪水槽。Phung[2]基于VOF法建立兩相流數(shù)值模型，并將其應(yīng)用于波浪與潛堤的相互作用以及波浪在斜坡上破碎的數(shù)值模擬。Lin[3]通過(guò)求解雷諾時(shí)均方程，采用k-ε模型模擬紊流作用，實(shí)現(xiàn)了孤立波和橢圓余弦波的數(shù)值模擬。王永學(xué)[4]采用VOF法捕捉自由表面，利用雷諾時(shí)均模型模擬波浪在直立墻前的破碎過(guò)程，得到直墻面上波壓力。齊鵬和王永學(xué)[5]基于VOF法建立三維數(shù)值波浪水槽，設(shè)置動(dòng)網(wǎng)格造波邊界，較好地模擬了規(guī)則波和結(jié)構(gòu)物的相互作用。李勝忠[6]基于FLUENT軟件建立二維數(shù)值波浪水槽，很好地完成了線性規(guī)則波、Stokes波、孤立波以及不規(guī)則波的模擬。

但無(wú)論是物理實(shí)驗(yàn)還是數(shù)值模擬，受到波浪水槽長(zhǎng)度和寬度影響，造波板產(chǎn)生的行進(jìn)波在遇到結(jié)構(gòu)物后會(huì)發(fā)生一次反射，一次反射波與原行進(jìn)波疊加后傳遞到造波板上會(huì)發(fā)生二次反射。因此在進(jìn)行一段時(shí)間的造波后，水槽內(nèi)波浪的多次反射、疊加將嚴(yán)重影響實(shí)驗(yàn)精度。Salter[7]等以搖板板上所受波壓力為反饋信息，設(shè)計(jì)遞歸式濾波器，對(duì)造波板的速度進(jìn)行調(diào)整，以解決二次反射的問(wèn)題。Bullock 和 Murton[8]以板上波高信號(hào)為反饋信息，不考慮造波板前的衰減駐波和波浪非線性的影響設(shè)計(jì)了遞歸式濾波器，得到造波板的位移，在物理水槽中實(shí)現(xiàn)了楔形造波機(jī)規(guī)則波和不規(guī)則波的主動(dòng)吸收式造波測(cè)試。國(guó)內(nèi)方面，陳漢寶[9]等根據(jù)勢(shì)函數(shù)計(jì)算造波機(jī)參數(shù)，探討了無(wú)反射造波機(jī)的實(shí)現(xiàn)思路。柳淑學(xué)[10-11]等提出代表頻率時(shí)域控制法用于不規(guī)則波的主動(dòng)吸收，實(shí)現(xiàn)了推板式造波機(jī)不規(guī)則波的主動(dòng)吸收式造波實(shí)驗(yàn)。顧民[12]等采用頻域?yàn)V波器方法設(shè)計(jì)完整具有主動(dòng)吸收功能的造波裝置，實(shí)現(xiàn)了搖板式造波機(jī)規(guī)則波和不規(guī)則波主動(dòng)吸收式造波實(shí)驗(yàn)。李宏偉[13]在時(shí)域控制法和頻域?yàn)V波器方法的基礎(chǔ)上，設(shè)計(jì)了基于S面控制的主動(dòng)吸收造波系統(tǒng)。楊洪齊[14]在時(shí)域和頻域兩個(gè)方向?qū)崿F(xiàn)了規(guī)則波和不規(guī)則波的主動(dòng)吸收，提出了帶有延時(shí)補(bǔ)償?shù)闹鲃?dòng)吸收控制方程，提高了造波質(zhì)量。

將濾波器方法應(yīng)用在主動(dòng)吸收式造波中，可以針對(duì)不同頻率的波浪做出不同的反應(yīng)效果，對(duì)不規(guī)則波有著更好的吸收效果。但同時(shí)由于濾波器方法對(duì)頻率的敏感特性，很容易受到橫波、不規(guī)則雜波等問(wèn)題影響，經(jīng)常會(huì)出現(xiàn)造波板異常抖動(dòng)、造波質(zhì)量較差等問(wèn)題。通過(guò)數(shù)值模擬方法可以更好地驗(yàn)證設(shè)計(jì)濾波器主動(dòng)吸收效果，但目前主動(dòng)吸收數(shù)值波浪水槽[15-16]大多采用時(shí)域控制主動(dòng)吸收算法，濾波器方法在數(shù)值模型中很少涉及。本文基于開(kāi)源軟件OpenFOAM，利用動(dòng)網(wǎng)格技術(shù)模擬推板式造波，自定義濾波器主動(dòng)吸收式造波邊界，建立粘性水槽數(shù)值模型驗(yàn)證濾波器方法的準(zhǔn)確性和可靠性，為物理模型實(shí)驗(yàn)提供參考。

1 數(shù)值水槽基本方法及理論

1.1 黏性波浪數(shù)值水槽的控制方程

采用VOF方法追蹤自由表面，其原理是定義項(xiàng)體積分?jǐn)?shù)α，當(dāng)α=0時(shí)表示該網(wǎng)格內(nèi)均為空氣，當(dāng)α=1時(shí)表示該網(wǎng)格內(nèi)均為水，當(dāng)0<α<1 時(shí)表示為空氣和水的交界面。二相流的平均密度和動(dòng)力粘性系數(shù)為

ρ=αρ1+(1-α)ρ2

(1)

μ=αμ1+(1-α)μ2

(2)

故黏性數(shù)值水槽的控制方程包括連續(xù)性方程、包括項(xiàng)體積分?jǐn)?shù)的動(dòng)量方程和體積分?jǐn)?shù)函數(shù)。

連續(xù)性方程

·U=0

(3)

動(dòng)量方程

(4)

體積分?jǐn)?shù)函數(shù)

(5)

1.2 推板造波理論

(6)

根據(jù)線性疊加法，不規(guī)則波波面可由有限個(gè)規(guī)則波波面線性疊加而成

(7)

N為規(guī)則波的數(shù)量，φi為隨機(jī)初相位，位于[0，2π]內(nèi)，Ai為第i個(gè)規(guī)則波的振幅，其表達(dá)式為

(8)

式中：ωi為第i個(gè)組成波的代表頻率；△ωi為頻率間隔；S(ωi)為代表頻率為ωi下的譜能密度。

根據(jù)線性造波理論，不規(guī)則造波時(shí)推板式造波板的位移為

(9)

1.3 頻域?yàn)V波器主動(dòng)吸收理論

造波板運(yùn)動(dòng)產(chǎn)生的波浪傳播到建筑物上發(fā)生一次反射，一次反射波與原行進(jìn)波疊加后傳遞到造波板上會(huì)發(fā)生二次反射，所以在進(jìn)行一段時(shí)間的造波后，水槽內(nèi)的波浪經(jīng)過(guò)多次反射、疊加后對(duì)試驗(yàn)精度產(chǎn)生較大影響。主動(dòng)吸收式造波是在造波板原有位移上添加一個(gè)修正位移來(lái)消除二次反射波，在水槽內(nèi)得到穩(wěn)定波場(chǎng)。造波板的位移可以用下式表示

x=xp+xa

(10)

造波板上波面為

η0=ηp+ηa+ηps+ηas+ηr+ηrr

(11)

式中：η0為造波板上實(shí)測(cè)的總波高；ηp和ηa為由xp和xa所產(chǎn)生的行進(jìn)波；ηps和ηas為xp和xa所產(chǎn)生的非傳播模態(tài)項(xiàng)；ηr為水槽中的反射波；ηrr為在造波板上產(chǎn)生的二次反射波。在滿足主動(dòng)吸收的條件下，有

ηa+ηrr=0

(12)

由此聯(lián)立方程推導(dǎo)出頻域傳遞函數(shù)為

(13)

式中：C0為水動(dòng)力傳遞系數(shù)；ECn為推板前駐波影響因子，分別與波數(shù)k0、虛波數(shù)kn和水深h有關(guān)。

得到的F即為主動(dòng)吸收式造波系統(tǒng)的主動(dòng)吸收傳遞函數(shù)，故濾波器法主動(dòng)吸收系統(tǒng)可以認(rèn)為是一個(gè)輸入為2ηp-η0，輸出為x，頻率響應(yīng)特性為F的單輸入、單輸出的線性時(shí)不變系統(tǒng)。

2 濾波器設(shè)計(jì)

對(duì)于一個(gè)典型的單輸入、單輸出的離散線性時(shí)不變系統(tǒng)，可用線性常系數(shù)差分方程來(lái)表示

(14)

將上式的N階線性常系數(shù)差分方程變換形式，取a0=1得

(15)

已知一個(gè)離散時(shí)間系統(tǒng)可用式(15)的差分方程表示，則在z域中系統(tǒng)函數(shù)可表示為

(16)

在主動(dòng)吸收系統(tǒng)中，輸入信號(hào)為目標(biāo)波高與實(shí)測(cè)波高的差值，是一個(gè)有界量；輸出信號(hào)為造波板位移，受到造波板幅值限制也是有界的。從定義上說(shuō)，如果輸入有界，則輸出必定有界的系統(tǒng)是穩(wěn)定的，因此頻域?yàn)V波器的設(shè)計(jì)也應(yīng)滿足穩(wěn)定性要求。但由于IIR濾波器系統(tǒng)函數(shù)存在分母項(xiàng)，故存在極點(diǎn)，有不穩(wěn)定的可能。從定義上看，系統(tǒng)穩(wěn)定的充要條件為系統(tǒng)的單位脈沖絕對(duì)可加和，即對(duì)于一個(gè)穩(wěn)定系統(tǒng)，存在正數(shù)B使得

式中：h(n)為系統(tǒng)的脈沖響應(yīng)。對(duì)于存在極點(diǎn)的濾波器系統(tǒng)，一般常用單位圓方法對(duì)系統(tǒng)穩(wěn)定性進(jìn)行判斷。當(dāng)極點(diǎn)落在單位圓外時(shí)，表示濾波器當(dāng)前輸出需要放大上一時(shí)刻的輸出，當(dāng)n→∞，系統(tǒng)的脈沖響應(yīng)趨于無(wú)窮大，不滿足絕對(duì)可加和條件。當(dāng)極點(diǎn)落在單位圓上時(shí)，當(dāng)n→∞，系統(tǒng)的脈沖響應(yīng)趨于常數(shù)或等幅值震蕩，不滿足絕對(duì)可加和條件。當(dāng)極點(diǎn)落在單位圓內(nèi)時(shí)，表示濾波器當(dāng)前輸出需要縮小上一時(shí)刻的輸出，當(dāng)n→∞，系統(tǒng)的脈沖響應(yīng)趨于零，滿足絕對(duì)可加和條件。對(duì)于有多個(gè)極點(diǎn)的系統(tǒng)，可根據(jù)代數(shù)學(xué)知識(shí)進(jìn)行分解，將其看做多個(gè)單極點(diǎn)并聯(lián)而成的系統(tǒng)。在并聯(lián)系統(tǒng)中，只要有一個(gè)系統(tǒng)不穩(wěn)定，整個(gè)系統(tǒng)就是不穩(wěn)定的。因此系統(tǒng)穩(wěn)定的充要條件為濾波器的極點(diǎn)均落在單位圓內(nèi)。

濾波器采樣頻率定義了濾波器每秒從輸入信號(hào)中提取并組成離散信號(hào)的采樣個(gè)數(shù)，所以濾波器的采樣頻率應(yīng)小于或等于輸入信號(hào)——波高傳感器的采樣頻率。本文以天津大學(xué)創(chuàng)新研究院的波流水槽作為工程背景，其采用波高儀的采樣頻率為1 000 Hz，適當(dāng)減小采樣頻率在實(shí)驗(yàn)室造波實(shí)驗(yàn)中會(huì)給上位機(jī)更多的時(shí)間對(duì)信號(hào)進(jìn)行處理，保證在每個(gè)時(shí)間間隔內(nèi)給出造波板的位移信號(hào)；在數(shù)學(xué)模型實(shí)驗(yàn)中可以適當(dāng)增大時(shí)間步長(zhǎng)、提高計(jì)算速度，因此本文濾波器設(shè)計(jì)頻率選為250 Hz。

本文通過(guò)求解色散方程得到不同水深下，C0、ECn與頻率的關(guān)系，進(jìn)而得到理論傳遞函數(shù)。采用基于線性最小二乘法(OLS)，選擇帶控制量的自回歸模型(ARX)作為擬合模型對(duì)理論傳遞函數(shù)進(jìn)行模擬。理論上，濾波器階數(shù)越高，擬合精度越高，但同時(shí)系統(tǒng)不穩(wěn)定、計(jì)算量增大、延遲問(wèn)題也會(huì)更加明顯。并且由于理論傳遞函數(shù)在低頻處幅值響應(yīng)十分大，如果對(duì)過(guò)低頻率的波浪進(jìn)行主動(dòng)吸收會(huì)導(dǎo)致造波板位移超過(guò)造波機(jī)限制，無(wú)法進(jìn)行正常造波。而在本文工程背景下，不需要對(duì)0.2 Hz以下規(guī)則波進(jìn)行主動(dòng)吸收，0.2 Hz以下波浪在不規(guī)則波的能量占比也較小，故根據(jù)實(shí)際需要，設(shè)計(jì)二階濾波器在低頻處對(duì)幅值進(jìn)行抑制，在常用實(shí)驗(yàn)造波頻率0.3～2.5 Hz范圍內(nèi)進(jìn)行擬合。

本文通過(guò)ARX方法得到在0.5 m水深下主動(dòng)吸收濾波器參數(shù)后首先對(duì)濾波器進(jìn)行穩(wěn)定性檢測(cè)，要保證濾波器所有的極值均落在單位圓內(nèi)。如圖1所示，其中×代表極點(diǎn)，○代表零點(diǎn)，兩極點(diǎn)坐標(biāo)分別為(0.998 4,0)和(0.955 4,0)，所有的極點(diǎn)均落在單位圓以內(nèi)，系統(tǒng)穩(wěn)定。其次要對(duì)濾波器使用要求進(jìn)行檢測(cè)，如圖2所示，濾波器在低頻處對(duì)幅值進(jìn)行抑制，在常用造波頻率范圍內(nèi)擬合程度良好，滿足實(shí)驗(yàn)要求。

圖1 濾波器零極點(diǎn)分布圖

3 自定義造波邊界

根據(jù)OpenFOAM自帶的網(wǎng)格運(yùn)動(dòng)邊界基類fixedValuePointPatchFields開(kāi)發(fā)新的濾波器主動(dòng)吸收式造波邊界類。此邊界包括私有數(shù)據(jù)時(shí)間序列、波面序列、濾波器階數(shù)和濾波器系數(shù)等。根據(jù)字典文件讀取關(guān)鍵字來(lái)得到以上私有數(shù)據(jù)。通過(guò)字典文件讀取時(shí)間-理論波面序列，得到當(dāng)前時(shí)刻的理論波面。通過(guò)VOF方法得到邊界上的水深，減去初始水深得到造波板上波面。根據(jù)濾波器階數(shù)和參數(shù)設(shè)計(jì)濾波器，以理論波面和實(shí)際波面為輸入量，通過(guò)儲(chǔ)存之前的位移及波面，根據(jù)IIR濾波器原理，利用差分方程計(jì)算得出當(dāng)前造波板位移。

4 數(shù)值水槽的建立及驗(yàn)證

4.1 數(shù)值水槽的設(shè)置

本文以天津大學(xué)創(chuàng)新研究院的波流水槽作為工程背景，如圖3所示，根據(jù)實(shí)際水槽條件建立一個(gè)長(zhǎng)9 m、高0.8 m、水深為0.5 m的二維矩形水槽模型。根據(jù)查晶晶[18]推薦,網(wǎng)格長(zhǎng)寬比不應(yīng)該超過(guò)20∶1，垂向波高內(nèi)網(wǎng)格數(shù)不應(yīng)少于10個(gè)，縱向波長(zhǎng)范圍內(nèi)網(wǎng)格不應(yīng)少于40個(gè)，故在X方向布置500個(gè)網(wǎng)格，在Z方向布置300個(gè)網(wǎng)格。左側(cè)為造波邊界，右端為全反射墻體邊界。根據(jù)濾波器的采樣頻率，時(shí)間步長(zhǎng)選為0.004 s。據(jù)水槽左端3.17 m、3.47 m、3.83 m處分別設(shè)置監(jiān)測(cè)點(diǎn)，得到該位置上的波面-時(shí)間曲線。在相同的波浪要素條件下，分別進(jìn)行主動(dòng)吸收造波和普通推板造波數(shù)值對(duì)比實(shí)驗(yàn)，對(duì)結(jié)果進(jìn)行分析處理，驗(yàn)證濾波器法主動(dòng)吸收的正確性。

圖3 數(shù)值波浪水槽示意圖

4.2 規(guī)則波的驗(yàn)證與對(duì)比

規(guī)則波數(shù)值模擬時(shí)長(zhǎng)定為60 s，目標(biāo)波形為余弦波。分析對(duì)比1號(hào)測(cè)點(diǎn)上主動(dòng)吸收造波和普通推板造波波面曲線。從圖4中可以看出主動(dòng)吸收式造波時(shí)，入射波在到達(dá)水槽右端后形成一次反射波，入射波和一次反射波疊加后傳遞到造波板時(shí)，造波板有效吸收了二次反射波，此時(shí)水槽內(nèi)只有入射波與一次反射波，二者疊加形成穩(wěn)定駐波。而采用普通推板式進(jìn)行造波，波浪發(fā)生多次反射，無(wú)法形成穩(wěn)定波場(chǎng)。證明了頻域?yàn)V波器方法對(duì)規(guī)則波的主動(dòng)吸收效果十分顯著。

4-a T=1.0 s,H=0.04 m 4-b T=1.5 s,H=0.04 m

4.3 不規(guī)則波的驗(yàn)證與對(duì)比

不規(guī)則波主動(dòng)吸收造波模擬時(shí)間為120 s，靶譜選為JONSWAP譜，根據(jù)1、2、3號(hào)測(cè)點(diǎn)上波面情況根據(jù)三點(diǎn)法進(jìn)行入反射分離得到入射波波面曲線。與普通造波時(shí)入射波波面對(duì)比如圖5所示。從圖中可以看出開(kāi)啟主動(dòng)吸收時(shí)波峰波谷較普通造波時(shí)小，即能量較正常造波時(shí)小。初步說(shuō)明頻域?yàn)V波器法主動(dòng)吸收能夠?qū)Χ畏瓷洳ㄟM(jìn)行有效吸收，波場(chǎng)內(nèi)波浪不會(huì)發(fā)生多次反射，能量不會(huì)反復(fù)疊加。

5-a T=1.0 s,H=0.04 m 5-b T=1.5 s,H=0.04 m

為進(jìn)一步說(shuō)明濾波器法主動(dòng)吸收式造波對(duì)二次反射波的吸收效果，分別對(duì)主動(dòng)吸收造波和普通造波入射波面進(jìn)行頻譜分析，得到的頻譜和目標(biāo)靶譜對(duì)比如圖6所示?？梢钥闯鲩_(kāi)啟主動(dòng)吸收時(shí)頻譜曲線與靶譜擬合較好，普通造波時(shí)譜峰密度遠(yuǎn)大于靶譜，說(shuō)明波浪多次反射疊加后能量較大、波場(chǎng)混亂。

6-a T=1.0 s,H=0.04 m 6-b T=1.5 s,H=0.04 m

為定量分析濾波器法對(duì)不規(guī)則波的吸收情況，參考王先濤[19]的分析方法定義主動(dòng)吸收率，根據(jù)主動(dòng)吸收造波前后頻譜曲線面積對(duì)比來(lái)判斷吸收效果。主動(dòng)吸收率ra計(jì)算公式如下

(17)

式中：Sp為正常造波頻譜曲線；Sa為主動(dòng)吸收式造波頻譜曲線；Sn為理論頻譜曲線。

同時(shí)分析對(duì)比特征周期、譜峰波高、譜峰密度來(lái)驗(yàn)證不規(guī)則波模擬精度。表1為濾波器法主動(dòng)吸收造波入射頻譜曲線的特征參數(shù)與靶譜特征參數(shù)的統(tǒng)計(jì)結(jié)果。不同周期下的特征周期T1/3、譜峰波高Hp和譜峰密度Sp誤差均小于3%。在1.5 s周期下，可能普通推板造波時(shí)能量較小導(dǎo)致吸收率與其余各組相比較低，為91.15%。其余各個(gè)周期下的吸收率均超過(guò)95%。綜合考慮各類波浪要素和吸收率，進(jìn)一步說(shuō)明濾波器法對(duì)不規(guī)則波有著良好的吸收效果。

表1 不規(guī)則波波浪要素統(tǒng)計(jì)結(jié)果

5 結(jié)論

本文基于OpenFOAM 不可壓縮二相流求解器interDyMFoam，開(kāi)發(fā)了仿物理造波數(shù)值水槽，實(shí)現(xiàn)了濾波器法主動(dòng)吸收式造波。通過(guò)繼承動(dòng)邊界類fixedValuePointPatchField開(kāi)發(fā)了濾波器法主動(dòng)吸收造波邊界。結(jié)合工程背景需要，確定了針對(duì)目標(biāo)水深下的IIR濾波器參數(shù)。

在不同周期規(guī)則波的造波實(shí)驗(yàn)中，濾波器法有效吸收了二次反射波，能夠得到穩(wěn)定波場(chǎng)。在不同周期不規(guī)則的造波實(shí)驗(yàn)中，波浪要素與靶譜目標(biāo)參數(shù)誤差均小于3%，濾波器法對(duì)二次反射波的吸收率均在90%以上。對(duì)比不同周期下規(guī)則波和不規(guī)則波造波效果，可以看出濾波器方法對(duì)規(guī)則波和不規(guī)則波的二次反射都有著良好的吸收效果。為后續(xù)濾波器方法在物理模型實(shí)驗(yàn)的研究提供了重要基礎(chǔ)依據(jù)。同時(shí)濾波器方法將主動(dòng)吸收問(wèn)題從時(shí)域映射到了頻域內(nèi)，可以通過(guò)對(duì)理論函數(shù)曲線擬合的不斷優(yōu)化，為解決實(shí)際造波問(wèn)題提供了更多可行的方向，也為解決三維水池主動(dòng)吸收造波問(wèn)題做好鋪墊。