緩坡方程的有限分析數(shù)值模式*

黃武平,張庭榮

(廣東省水利水電科學(xué)研究院,廣東省水動力學(xué)應(yīng)用研究重點實驗室，廣東 廣州 510635)

?

緩坡方程的有限分析數(shù)值模式*

黃武平,張庭榮

(廣東省水利水電科學(xué)研究院,廣東省水動力學(xué)應(yīng)用研究重點實驗室，廣東 廣州 510635)

基于時間關(guān)聯(lián)的雙曲型緩坡方程,采用九點有限分析法,建立近岸波浪傳播變形的數(shù)值模式,在對時間關(guān)聯(lián)的雙曲型緩坡方程數(shù)值離散時,空間導(dǎo)數(shù)采用九點有限分析法,時間導(dǎo)數(shù)仍采用有限差分法,再采用質(zhì)量集中方法得到本文的數(shù)值模式。為了驗證本模式的數(shù)值精度和有效性,分別對Berkhoff(1972)實驗和突堤地形進行數(shù)值模擬,結(jié)果表明本模式的模擬值與物理模型試驗值、解析解及有限差分模式的模擬值吻合良好,表明本模式可以對近岸波浪傳播過程中折射、繞射、反射、淺水變形、波浪破碎和弱非線性等因素的影響進行較好的預(yù)測。

緩坡方程;時間關(guān)聯(lián);有限分析法;質(zhì)量集中法;數(shù)值模擬

有限分析法(FAM)是美國IOWA大學(xué)陳景仁教授于20世紀80年代初提出的一種數(shù)值計算方法[1-3]。FAM是有限元基礎(chǔ)上的一種改進方法,與有限元方法和有限差分方法相似,FAM也需要把計算區(qū)域劃分為單元子域。他的基本思想是：在離散單元內(nèi)求控制方程的局部分析解,如方程為非線性的,求分析解時,在子域內(nèi)將其線性化。用有限分析法解出的數(shù)值解是九點格式的,具有精度高,明顯的自動迎風性質(zhì),計算穩(wěn)定性好等優(yōu)點。

1 數(shù)值求解

1.1 有限分析法求解

所有的數(shù)值計算方法,都是把連續(xù)解變成離散解,從而把偏微分方程變成代數(shù)方程,利用電子計算機求解。各種方法的差異在于如何建立所求點(中心點)和周圍點之間的關(guān)系。二維定常納維爾-斯托克斯方程為：

LU=Re·(UUx+VUy)-(Uxx+Uyy)=G

(1)

其中L表示解算子,G表示重力或壓力,表示雷諾數(shù),U、V分別表示沿x和y方向的流速值。

有限分析法基本思想是：在離散單元內(nèi)的解U不再用插值函數(shù)式來表達,而是方程局部線性化后的分析解。以方程式(1)為例,在單元內(nèi)局部線性化后的方程寫為：

(2)

圖1 有限分析法

1.2 緩坡方程有限分析數(shù)值模式的建立

對于時間關(guān)聯(lián)型的雙曲型緩坡方程,把時間變量看成常量時,對于空間變量來說,時間關(guān)聯(lián)型的雙曲型緩坡方程變成泊松方程形式,再運用有限分析法對緩坡方程求數(shù)值解,得到本文的數(shù)值模式。

Berkhoff (1972)[4]在線性理論中引入一個表示地形緩變的小參數(shù),并利用攝動法和沿水深積分的方法導(dǎo)出了一個二階偏微分方程式,稱之為緩坡方程。

▽·(CCg▽Ψ)+k2CCgΨ=0

(3)

其中 速度勢函數(shù)φ(x,y,t)和Ψ(x,y)的關(guān)系為：φ(x,y,t)=Re{Ψ(x,y)e-iωt},C和Cg分別為行進波的局部相速度和局部群速度。

(4)

▽2ζ=-ζ

(5)

如果網(wǎng)格尺寸均勻,即x,y兩個方向的網(wǎng)格各自均勻分布,運用上述有限分析法求解可以得出方程(5)的解：

(6)

為了得到顯式數(shù)值格式,對特解部分采用質(zhì)量集中法進行簡化,得到如下形式：

(7)

將物理量(ζ,k,C,Cg和W)布置在網(wǎng)格節(jié)點上,對空間步長導(dǎo)數(shù)Δx,Δy和時間步長導(dǎo)數(shù)Δt進行離散,空間導(dǎo)數(shù)離散采用有限分析法,時間導(dǎo)數(shù)離散仍采用有限差分,方程(7)可離散為:

(8)

2 數(shù)值模式的檢驗與比較

2.1 Berkhoff實驗

Berkhoff et al在斜坡和橢圓型淺灘組合地形上進行了波浪傳播的物理模型試驗,試驗地形如圖2圖3所示,淺灘的中心位于(11.0 m,10.0 m),入射波周期T為1.0 s,波高H0是0.05 m,1#～8#斷面分別位于x=12,14,16,18,20和y=8,10,12 m。這是個經(jīng)典算例,用以檢驗?zāi)Ｊ綄Σɡ苏凵洹⒗@射和非線性的有效性。

模式的計算參數(shù)如下：計算域為22 m×20 m,左右邊界為全反射邊界,下游邊界為開邊界。模型的空間步長Δx=Δy=0.1 m,時間步長Δt=T/40,計算時間為60T。由于計算域相對較小,所以不考慮海底摩擦的影響,即fb=0。

圖2 Berkhoff地形水深等值線

圖3 試驗地形剖面示意

圖4是數(shù)值模擬的A、B、C、D、E、F六點波面隨時間變化的過程線；圖5是數(shù)值模擬的斷面在t=60T時刻的瞬時波面圖；圖6是1#～8#斷面相對波高的數(shù)值解與物理模型實驗值的比較；圖4-11是數(shù)值模擬的波面平面分布圖。圖4和5表明,模式的數(shù)值結(jié)果穩(wěn)定、收斂；圖6表明,兩種模式的數(shù)值解與物理模型實驗值吻合較好,非線性擴展模型的數(shù)值解比線性擴展模型數(shù)值解更接近實驗值,尤其在淺灘后非線性作用較強區(qū)(見圖6中的3#～8#),而且九點有限分析格式模式比五點有限差分格式模式模擬的效果略好(見6中的6#)；圖7表明,在橢圓型淺灘后,由于波浪的折射作用形成一了個波能輻聚的主峰,同時由于繞射作用使波能沿波峰線側(cè)向傳遞,這樣就在橢圓型淺灘后形成了一系列的波能輻聚和輻散區(qū)域。表明本模式對波浪的折射、繞射和非線性有較高的有效性。

圖4 數(shù)值模擬的 A、B、C、D、E、F六點波面隨時間變化的過程線

圖5 數(shù)值模擬7#斷面在t=60T時刻的瞬時波面示意

圖6 九點有限分析法和五點有限差分法相對波高比的數(shù)值解和物理模型實驗值的比較

圖7 數(shù)值模擬的波面平面分布示意

2.2 突堤附近波浪場

為了檢驗本模式對斜向波的折射、繞射、反射、淺水變形、波浪破碎和非線性的實用性,本文對突堤附近波浪場進行了數(shù)值模擬,將計算結(jié)果與關(guān)孟儒模型試驗的結(jié)果進行比較[6]。物理模型實驗區(qū)域為6 m×10 m,地形為1∶50的斜坡,突堤為從靜止水時水邊線起長4 m的薄板形狀,周期為1.2 s、波高為2 cm的入射波以30°角度從左邊邊界入射。計算區(qū)域如圖8所示,計算域?qū)?0 m,長5.75 m(右邊界的水深為0.5 cm,其附近波浪已破碎),上下邊界及防波堤兩側(cè)均為全反射邊界。空間步長Δx=Δy=0.025 m,時間步長Δt=0.0012 s,底摩擦因數(shù)fb=0.01。

圖9是模式數(shù)值模擬的波高等值線圖,圖10是模式數(shù)值模擬的波面平面分布圖,圖11是不同斷面沿岸灘方向的波高值與實驗值的比較。從圖9～圖10可知,在斜向波影響下突堤前面出現(xiàn)干涉區(qū),突堤后面出現(xiàn)繞射區(qū),前進波的波高隨水深的變淺逐漸變大,最后由于波浪破碎,波能衰減,岸灘處波高減小,圖9中實心圓點為物理模型記錄的波浪破碎點,這與波高等值線相吻合；從圖11可知,非線性擴展模型比線性擴展模型模擬的效果更好,計算值更接近實驗值,突堤前反射區(qū)數(shù)值模擬計算波高值幾乎達到入射波高的2倍,比實驗值要偏高,尤其在突堤的頂端處更顯著,這是由于數(shù)值計算時突堤頂端發(fā)生全反射,突堤后面的繞射區(qū)的數(shù)值模擬計算值跟實驗值比較吻合,離突堤比較近的斷面(y=4.8 m)處繞射比較弱,而離突堤較遠的斷面(y=3 m)處繞射比較強。表明本模型對斜向波的折射、繞射、反射、淺水變形、波浪破碎和弱非線性效應(yīng)有比較好的實用性。

圖8 計算區(qū)域示意

圖9 本模式數(shù)值模擬的波高等值線示意

圖10 本模式數(shù)值模擬的波面平面分布示意

圖11 不同斷面波高分布示意

3 結(jié)語

本文以趙紅軍提出的時間關(guān)聯(lián)的雙曲型擴展緩坡方程為控制方程,采用九點有限分析法和質(zhì)量集中法建立了近岸波浪場的有限分析數(shù)值模式。運用本模式數(shù)值模擬近岸波浪場的波浪傳播,通過對Berkhoff地形實驗和突堤進行數(shù)值模擬檢驗了本模式對波浪的折射、繞射、反射、淺水變形、波浪破碎和非線性有很好的實用性。

[1] Chen C J，Chen H C. Finite analytic numerical method for unsteady two-dimendional navier-stokes equations[J]. Journal of Computational Physis，1984(53)：209-226.

[2] Chen C J. Finite analytic numerical solution of heat transfer in two-dimensional cavity flow [J]. Numerical Heat Transfer，1981(4)：179-197.

[3] 陳景仁. 流體力學(xué)及傳熱學(xué)[M].北京：國防工業(yè)出版社，1984.

[4] Berkhoff J. C. W.． Computation of combined refraction-diffraction[C]∥Proceedings of the 13th International Conference on Coastal Engineering，ASCE 1972：745-747．

[5] 趙紅軍，宋志堯，徐福敏，等. 一種擴展型緩坡方程的時域計算模式[J].水動力學(xué)研究與進展A輯，2009,24(4)：503-511.

[6] 關(guān)孟儒，牛克源，湯乃銘，等. 海洋動力學(xué)導(dǎo)論[M]. 南京：河海大學(xué)出版社，1995.

(本文責任編輯 王瑞蘭)

Finite Analytical Numerical Mode for Mild Slope Equation

HUANG Wuping, ZHANG Tingrong

(Guangdong research institute of water resources and hydropower, Guanghou 510635, China)

A numerical model for wave propagation is proposed, which is accomplished by nine point finite-analytic methods based on the time-dependent hyperbolic mild-slope wave equation. The space’s differential is introduced by nine point finite-analytic method and the time’s differential is introduced by finite-difference method when the time-dependent hyperbolic mild-slope wave equation is dispersed, and the numerical model has been got by the lumped mass method. The Berkhoff(1972) experiment and a breakwater stand out are simulated for validating the precision and validity of my numerical model. The calculated results are in agreement with the experimental value, theoretical and finite-difference model’s solution, that indicate that my model is capable of giving good predictions of wave refraction, diffraction, shoaling, breaking energy dissipation and weakly nonlinearity in the near shore zone.

mild-slope equation, time-dependent, finite-analytic method, the lumped mass method, numerical simulations

2016-07-08;

2016-07-15

廣東省水利科技創(chuàng)新項目(2012-06)。

黃武平(1985),男,碩士,工程師,從事水動力學(xué)研究。

TV139.2

A

1008-0112(2016)08-0001-06