耦合算法在幕墻式消浪結構性能研究中的應用

1、耦合算法在幕墻式消浪結構性能研究中的應用             摘要作者針對近岸波浪與結構物相互作用問題提出了一種耦合數(shù)值方法，即用時均化的二維雷諾平均的Navier stokes方程-流體體積法模型表達內域流動，用一維Boussinesq方程表達外域流動，通過速度、壓力和波面匹配邊界條件實現(xiàn)兩種數(shù)值模型的同步求解。耦合模型中的二維子模型能夠較好地表達結構物附近流動的細部特征，包括漩渦結構；一維子模型的計算效率很高，可通過延長其計算域以達到有效地避免二次反射波的影響。所建立

2、的耦合數(shù)值模型被證實可應用于幕墻式消浪結構防波性能的研究。 關鍵詞幕墻式消浪結構 耦合數(shù)值方法 VOF方法 Boussinesq方程 以往的研究成果表明，迎浪面開孔的沉箱直立堤可以有效地減小反射波，但消浪室的寬度(即開孔前墻和不透水后墻之間的寬度)一般應達到當?shù)夭ㄩL的四分之一1。如果入射波為涌浪或者其他類型的長波，這意味著理想的消浪室寬度在實際工程上可能無法實現(xiàn)。最近，日本學者提出了一種能有效消減直立堤前反射波的新型結構幕墻式消浪結構(curtain-walled dissipater)，其斷面如圖1所示。設在直立墻前的垂直屏障稱為幕墻。幕墻至直立墻的距離B為消浪室的寬度。圖1 幕墻式消浪結構

3、幕墻吃水深度用c表示。由于入射波引起的消浪室內水體振蕩運動和幕墻下面的渦旋運動使得波能大量耗散，從而實現(xiàn)消減反射波浪的目的。這種新型消浪結構的主要優(yōu)點是能夠有效地減小消浪室的設計寬度。數(shù)值模擬是揭示幕墻式消浪結構水力學性能和消浪機理的有效手段。由于在幕墻和直立墻處產(chǎn)生的反射波在造波邊界處可能形成二次反射，通常的方法需要在二次反射波傳播到結構物之前停止計算。這意味著計算域的長度必須足夠大。然而，在一個很長的立面二維計算域上全部采用粘流波浪數(shù)值模型做精細模擬一方面計算工作量很大，另一方面必要性也不充分。為此，本文提出了一個耦合求解策略，即將二維RANS-VOF模型與一維Boussinesq方程模型

4、耦合起來解決問題。在耦合模型中，一維子模型由于其計算效率很高，可以考慮足夠長的計算域；二維子模型則能夠較好地反映流場的細部，包括粘性對流動的影響。1 耦合模型的原理如圖2所示，耦合模型是將整個計算域劃分為1和2兩個子域。這兩個子域通過一條公共的重疊帶銜接起來。1為包含幕墻和直立墻的近場，流動以二維紊流運動方程，即Reynolds方程(RANS)為控制方程，采用標準k-紊流模型封閉，并在近壁區(qū)應用壁面函數(shù)理論；自由水面的描述采用定義流體體積函數(shù)的方法。圖2 耦合模型區(qū)域劃分 2域內流動的控制方程采用色散性改進的Boussinesq方程的一維形式，經(jīng)差分離散后得到系數(shù)矩陣為三對角矩陣形式的代數(shù)方程

5、組，采用追趕法快速求解。在耦合模型中，兩個子模型RANS和Boussinesq各自獨立求解，耦合的實現(xiàn)體現(xiàn)在重疊帶上流動信息的匹配。為了便于耦合處理，Boussinesq方程和RANS方程均采用交錯網(wǎng)格進行差分離散。其中，Boussinesq方程的離散參考了Madsen和Sorensen所用的格式。RANS動量方程中時間項的離散格式為向前差分，粘性項的離散格式為二階中心差分。為消除數(shù)值粘性的影響，動量方程中對流項的離散格式采用了三階迎風差分格式。差分方程的求解采用了SOLA-VOF方法。其基本思想是：首先用前一時刻的流場計算結果代入動量方程的顯式差分格式，求出當前時刻流場的近似值；再通過對壓力

6、廚行迭代修正，使得連續(xù)方程在一定的精度條件下得以滿足，對表面單元要求滿足自由表面的動力學邊界條件，即通過線性插值確定表面單元中心處的壓力值；在完成壓力迭代后，再對速度進行校正，然后用校正后的速度值代入k-方程相應的差分格式求解紊動動能和紊動動能耗散率；最后，應用施主與受主單元模型計算當前時刻的流體體積函數(shù)，確定流體自由表面的位置。由于動量方程、紊動動能方程和紊動動能耗散率方程對近壁區(qū)網(wǎng)格細密程度的要求不同，耗散率方程的要求最嚴，動量方程和動能方程的要求基本一致，為了既保證解的精度而又不致使網(wǎng)格劃分太密，本文在固壁區(qū)附近采用了壁面函數(shù)方法進行處理。即在壁面附近引入以下關系其中0是常數(shù)，與壁面粗糙

7、度有關，本文取0=0.0005；L是特征長度，計算中取為近壁區(qū)網(wǎng)格中心到壁面的距離。 合理地設置匹配邊界條件，使得內域流動和外域流動在匹配邊界處光滑而連續(xù)地過渡，是保證子模型耦合的關鍵。本文的做法如圖3所示。匹配邊界B是2域的出流邊界，Boussinesq模型執(zhí)行每一時間步的計算之前需要預知該邊界上的速度和波面值，由于邊界B同時又在1域的計算節(jié)點上，于是B邊界的速度條件B可利用1域得到的流場信息表達如下，圖3 匹配邊界附近的差分網(wǎng)格單元 (1)式中：u為水平速度，F(xiàn)為流體體積函數(shù)，y為垂向的網(wǎng)格步長，j為垂向網(wǎng)格節(jié)點編號，Jmax表示垂向網(wǎng)格節(jié)點的最大編號。給B邊界的波面賦予匹配條件時數(shù)值試驗

8、表明，需要利用連續(xù)方程反映的水位流量關系給出匹配條件其效果好于直接給定水位過程條件。這是因為通過水位和流量的相互調整，計算域內的反射波可得以減弱。因而，在實際計算中B邊界的波面可表示為利用Boussinesq方程推導過程中引入的近似展開關系，在波面函數(shù)和(深度平均)速度導的結果，V邊界上水平速度u“,垂向速度v和壓力p的匹配條件可按下面表達式給出耦合模型同步求解過程可簡單概括如下。首先，考慮2域左邊界處的入射波條件，在2域內執(zhí)行Boussinesq模型，當2域中接近匹配邊界的節(jié)點上的水平速度值第一次達到10-3m/s量級時，開始在1域內執(zhí)行RANS-VOF模型。在某一時間步n,執(zhí)行Boussi

9、nesq模型所需要的匹配邊界B處的速度和波面條件按式(1)和式(2)給出；Boussinesq模型在當前時間步的計算完成后，隨即利用式(3)(5)計算出匹配邊界V處的波面、速度和壓力邊界條件，并啟動RANS-VOF模型；RANS-VOF模型在當前時間步的計算完成后，即按照式(1)和(2)計算出匹配邊界B所在位置的速度和波面匹配條件。然后進入下一時間步的計算。         2 計算結果為了驗證計算結果，參考已知的模型試驗條件，取耦合模型為等水深(d=0.42m)，計算域全長18.36m.其中，子區(qū)域2和1的長度分

10、別為12.58m和5.78m.1域中幕墻厚度q=0.042m，幕墻下端為一45°尖角。幕墻吃水深度按相對吃水c/d=0.30.7計算了多個工況。幕墻和直立墻均作為全反射邊界處理。消浪室寬度B=0.29m取為固定值。通過改變入射波的周期，消浪室的相對寬度在B/L=0.080.16(L為波長)的范圍內變化。計算網(wǎng)格的劃分情況如下。2域按等步長(x=0.042m)劃分為300個單元格；1域的垂向也按等步長(z=0.021m)劃分，但其水平方向是按變步長劃分的，共劃分了150個單元格。1域中單元編號i=1120的各單元為xi=0.042m,i=121125的各單元為xi=0.035m，i=1

11、26130的各單元為xi=0.028m,i=131150的各單元為xi=0.021m.幕墻附近和消浪室區(qū)域的網(wǎng)格劃分得比較細的目的是為了更有效地分辨幕墻下面渦旋流動的特征。2域左邊界入射波的波高H=0.06m.應用合田方法測定反射率系數(shù)，在幕墻前大約2倍波長以外處設置了適當間隔的兩個浪高儀(如圖2所示)記錄波面變化。圖4為執(zhí)行耦合模型得到的瞬時(t=19.972s)波面(圖中點劃線為匹配邊界位置).圖5顯示了波浪反射系數(shù)Cr隨相對消浪室寬度B/L變化的關系。當幕墻吃水深度c/d=0.5一定時，Cr隨B/L變化，在B/L=0.12附近出現(xiàn)了最小值。本文的數(shù)值結果同模型試驗的結果十分接近，說明1D

12、/2D耦合模型是有足夠精度的。計算過程中觀察到了消浪室內水體做整體升降振蕩的現(xiàn)象。取消浪室平均水位變化的幅值(HT)與原入射波高(H)之比描述波浪在消浪室中被激勵的情況。圖6顯示了波浪在消浪室中水面波動被放大的比率隨消浪室相對寬度B/L變化的情況。可以看出，大約從B/L=0.12開始，消浪室中波浪的放大率隨B/L的增大而迅速減小的特點十分明顯。這說明，在給定幕墻吃水深度(c=0.12m)和消浪室寬度(B=0.29m)的情況下，波長(L)相對較大的波浪經(jīng)垂直幕墻透射到消浪室后引起的水體振蕩幅度更大。圖7顯示了幕墻反射系數(shù)Cr隨其相對吃水深度c/d變化的關系。當消浪室相對寬度B/L=0.12一定時

13、，Cr隨c/d變化，在c/d=0.5附近出現(xiàn)了最小值。H=0.06m，B/L=0.12，c/d=0.5圖4 整個計算域瞬時波面 d=0.42m,c=0.21m,H=0.06m圖5 Cr-B/L關系 d=0.42m,c=0.21m,H=0.06m圖6 消浪室中波浪放大率 B/L=0.12，H=0.06m圖7 Cr-c/d關系 對于幕墻反射系數(shù)為最小的情形(B/L=0.12，c/d=0.5)，圖8顯示了消浪室內的平均水位變化。從計算中可以看出，消浪室內水體的運動呈整體活塞式的運動。幕墻末端附近的復雜流態(tài)對自由水面并沒有大的影響。圖9顯示了幕墻附近渦旋流動的演化過程?？梢钥闯?，消浪室內常駐著一個逆時

14、針方向的漩渦，其尺度明顯受控于消浪室寬度。在消浪室水位由高向低變化的過程中，幕墻開口處外側有一個順時針方向的漩渦產(chǎn)生；隨著消浪室內水位的回升，該漩渦逐漸進入消浪室。 圖8 消浪室內平均水位變化(H=0.06m,B/L=0.12，c/d=0.5) 受消浪室內常駐漩渦的影響，進入消浪室內的漩渦有尺度減小的趨勢。消浪室水位再一次由高向低變化時，幕墻開口處外側將有新的順時針方向的漩渦產(chǎn)生。如此反復，渦漩的頻繁產(chǎn)生和移動導致了波浪能量的大量耗散，從而實現(xiàn)了消減反射波浪的目的。3 結語以上針對近岸波浪與結構物相互作用問題提出了一種耦合數(shù)值方法，即用二維RANS-VOF模型表達內域流動，一維Boussinesq方程模型表達外域流動，通過速度、壓力和波面匹配邊界條件實現(xiàn)兩個子模型的同步求解。耦合模型中二維子模型能夠較好地表達結構物附近流動的細部特征，包括漩渦結構；一維子模型的計算效率很高，使得大域問題的求解成為可能。所建立的耦合方法能有效地用于研究幕墻式消浪結構的水力性能。計算結果表明：(1)消浪室相對寬度B/L和幕墻相對吃水深度c/d都是影響垂直幕墻波浪反射系數(shù)Cr的重要參數(shù)。當c/d=0.5一定時，在B/L=0.12附近Cr為最小；當B/L=0.12一定時，在c/d=0.5附近Cr為最小。也即，當消浪室寬度接近0.12