基于濾波器方法主動吸收式數(shù)值水槽研究

劉 宣，白志剛，余海濤，高樹軍

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

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

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

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

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

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

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

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

(1)

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

(2)

故黏性數(shù)值水槽的控制方程包括連續(xù)性方程、包括項體積分數(shù)的動量方程和體積分數(shù)函數(shù)。

連續(xù)性方程

·U=0

(3)

動量方程

(4)

體積分數(shù)函數(shù)

(5)

1.2 推板造波理論

(6)

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

(7)

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

(8)

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

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

(9)

1.3 頻域濾波器主動吸收理論

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

x=xp+xa

(10)

造波板上波面為

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

(11)

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

ηa+ηrr=0

(12)

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

(13)

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

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

2 濾波器設(shè)計

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

(14)

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

(15)

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

(16)

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

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

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

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

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

圖1 濾波器零極點分布圖

3 自定義造波邊界

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

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

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

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

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

4.2 規(guī)則波的驗證與對比

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

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

4.3 不規(guī)則波的驗證與對比

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

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

為進一步說明濾波器法主動吸收式造波對二次反射波的吸收效果，分別對主動吸收造波和普通造波入射波面進行頻譜分析，得到的頻譜和目標靶譜對比如圖6所示。可以看出開啟主動吸收時頻譜曲線與靶譜擬合較好，普通造波時譜峰密度遠大于靶譜，說明波浪多次反射疊加后能量較大、波場混亂。

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

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

(17)

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

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

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

5 結(jié)論

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

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