曾維鴻，傅卓佳，湯卓超

（河海大學 力學與材料學院，南京 211100）

引 言

水波傳播廣泛存在于船舶和海洋工程領(lǐng)域.在水波傳播過程中，當水波到達近海岸地帶時，會引起海岸侵蝕，危及海岸建筑物的安全.在沿海地區(qū)，應(yīng)用水下防波堤保護近海岸建筑物，亟需準確預(yù)測水波與這些海上結(jié)構(gòu)物的相互作用[1].

在波浪傳播行為研究中，雖然實驗研究能夠直觀地模擬水波的傳播，但是由于其受到場地大小、耗時問題的困擾，使得實驗研究有著諸多限制.因此，數(shù)值模擬在波浪傳播力學行為研究中至關(guān)重要.經(jīng)典的數(shù)值方法包括有限差分法、有限單元法和邊界元法等，這些算法在數(shù)值水槽中得到了廣泛的應(yīng)用[2-6].例如： 王本龍等[2]采用有限差分法，初步模擬了波流的相互作用，為進一步研究波流相互作用奠定了基礎(chǔ).王大國等[3]基于有限單元法，對三維完全非線性數(shù)值波浪水槽進行了模擬.衛(wèi)志軍等[4]采用拓撲優(yōu)化技術(shù)研究了數(shù)值水槽中的晃蕩抑制.Koo 等[5]和Christou 等[6]使用邊界元法模擬了非線性水波與結(jié)構(gòu)物的相互作用行為.

然而，傳統(tǒng)網(wǎng)格類數(shù)值方法在求解此類移動邊界問題時往往面臨重復(fù)劃分網(wǎng)格導致的計算成本劇增難題.而基于無需網(wǎng)格劃分特點的無網(wǎng)格類配點方法在模擬波浪傳播問題時，比網(wǎng)格類方法具有固有和先天的優(yōu)勢.近些年來，有眾多無網(wǎng)格類方法被應(yīng)用到數(shù)值水槽計算中.例如，黃志濤等[7]采用光滑粒子流體動力學模擬了未滿載罐車液體晃蕩抑制；Senturk[8]采用徑向基函數(shù)法，對二維線性和非線性Stokes 波傳播進行了數(shù)值模擬；李珺璞等[9]采用奇異邊界法研究了水下障礙物對水波傳播的影響；在此基礎(chǔ)上，Zhang 等[10-11]和Huang等[12]采用廣義有限差分法（generalized finite difference method，GFDM），對數(shù)值水槽進行了模擬研究；Zhang 等[13]采用廣義有限差分法結(jié)合光滑粒子流體動力學，數(shù)值模擬了弱可壓縮黏性流的動力特性.

另一方面，一類半解析配點技術(shù)——邊界節(jié)點法[14]被提出，并成功用于聲波傳播和功能梯度材料熱傳導分析[15].由于該方法采用了預(yù)先滿足控制方程的非奇異半解析基函數(shù)，因此僅需少量離散節(jié)點即可得到高精度的計算結(jié)果.然而，類似于其他配點技術(shù)，傳統(tǒng)邊界節(jié)點法也會生成稠密矩陣，極大地限制了其在大規(guī)模問題中的應(yīng)用.近年來，Wang 等[16]和Xiong 等[17]先后將局部配點技術(shù)與傳統(tǒng)邊界節(jié)點法結(jié)合，建立了局部邊界節(jié)點法（localized boundary knot method，LBKM）.該方法是一種區(qū)域化離散方法，僅需要物理域內(nèi)的離散節(jié)點而無需任何網(wǎng)格.對于每個節(jié)點，可以首先通過節(jié)點之間的Euclid 距離確定具有簡單幾何的局部子域.然后，每個節(jié)點處的未知變量可以表示為該點對應(yīng)的局部子域內(nèi)節(jié)點處物理量的線性組合.將其代入控制方程和相應(yīng)的邊界條件后可以形成稀疏線性系統(tǒng)，其稀疏系統(tǒng)占用內(nèi)存小且計算時間短，這使得該方法更適合解決大規(guī)模問題.目前該方法已被成功用于聲學和對流擴散問題計算.

為了模擬波浪與防波堤的相互作用，防止近海岸建筑物發(fā)生侵蝕，本文首次應(yīng)用局部邊界節(jié)點法結(jié)合二階Runge-Kutta 法研究了波浪與水下防波堤之間的相互作用.第1 節(jié)介紹了數(shù)值波浪水槽的控制方程和邊界條件以及局部邊界節(jié)點法的求解方案；第2 節(jié)模擬了波浪與水下防波堤的相互作用；第3 節(jié)給出了結(jié)論.

1 數(shù)學模型

1.1 控制方程

本研究采用勢流原理進行波浪傳播的模擬，即假設(shè)液體無黏、無旋、不可壓縮，可以得到以下控制方程：

式 中，φ (x,z,t)為 速度勢，Ω為水槽計算域.

1.2 邊界條件

數(shù)值波浪水槽如圖1所示，在下邊界上，由于假設(shè)液體無法自由出入邊界，因此滿足不可穿透邊界條件：

圖1 數(shù)值波浪水槽的示意圖Fig.1 The schematic diagram of the numerical wave flume

式中，n為邊界上垂直向外的單位法向量.

數(shù)值波浪水槽主要是通過左側(cè)入射邊界條件施加入射波，從而達到波浪在水槽內(nèi)的傳播.由于初始水面為靜止狀態(tài)，為了逐步產(chǎn)生入射波從而增強水槽的穩(wěn)定性，所以引入斜坡函數(shù)Rm(t)，故入射邊界條件寫成

式中，U(z,t)為 入射波速度，Tm為調(diào)制時間，根據(jù)以往研究經(jīng)驗，本研究中取Tm=2T，T為入射波周期.

為了使到達右端邊界的波浪自由射出邊界，右端邊界采用輻射吸收邊界條件[11]：

式中，C為入射波的相速度，對于穩(wěn)定波而言C=λ/T，λ為入射波波長.

為了防止右側(cè)邊界反射波再次進入水槽，從而影響入射波的傳播，本文將采用海綿層和輻射吸收邊界進行波浪的逐步吸收，在半Lagrange 法的思想上，自由液面需滿足動力邊界條件和運動邊界條件[6]：

式中，η為自由表面高程函數(shù)，ν (x)為 阻尼系數(shù)，αs和 β 為調(diào)諧因子和長度因子，本文都取成常數(shù)1，b為水槽長度，ω為圓頻率.

1.3 局部邊界節(jié)點法

局部邊界節(jié)點法是一種新型的區(qū)域型無網(wǎng)格配點技術(shù)，僅需要在計算域內(nèi)布置離散點而無需任何網(wǎng)格.對于計算域中的離散點，都能夠找到包含本身和最鄰近的m個 點形 成的局部子域.為了統(tǒng)一，令x(i)=.對于局部子域中m+1個 相鄰節(jié)點，可以得到以下局部邊界節(jié)點法公式：

式中，rk j=‖xk?xj‖為 第k個點與第j個點之間的距離，(rk j)為微分算子的非奇異半解析基函數(shù).Laplace 算子的非奇異半解析基函數(shù)為平移不變調(diào)和函數(shù)：

其中，r=[r1,r2]，r1,r2分 別為中心點與其鄰近點在x方 向和z方 向上的距離；c為形狀參數(shù)，其取值與計算域特征長度有關(guān).

式（9）可以改寫成如下矩陣形式：

基于移動最小二乘原理的基本思想，需要在每個局部子域中定義以下殘差函數(shù)：

顯然使殘差Θ (u)取到最小值的 α(i)為偏微分方程的數(shù)值解.則可以對Θ (u)求 α(i)的變分，得到方程組

式中

因此，式（13）可以改寫為

把式（13）代入式（9）就可以得到x(i)處物理量的近似解：

同時，點x(i)處的一階和二階導數(shù)可以表示為

最后把式（18）和式（19）用于離散二階偏微分方程的控制方程和邊界條件，并求解相應(yīng)的代數(shù)方程組，從而近似得到偏微分方程的數(shù)值解.

1.4 自由邊界離散處理技術(shù)

由于式（6）和式（7）是與時間相關(guān)的非線性自由邊界條件，本文采用二階Runge-Kutta 法[18]分別離散自由液面邊界條件式（6）和式（7）中的時間項，其具體表達式如下：

式中，Δt為時間增量.當離散式（6）時，f與式（6）右端項相同，其中泛函V表示速度勢φ；類似地，當離散式（7）時，f與式（7）右端項相同，其中泛函V表示自由表面高程η.

2 數(shù)值波浪水槽特性研究

2.1 二階Stokes 波傳播

本節(jié)數(shù)值波浪水槽如圖1所示，長度b和 水深h分別為8 和1.若無其他特殊說明，本小節(jié)算例中采用的總節(jié)點數(shù)、時間間隔、鄰近點數(shù)固定為N=13 159, Δt=0．005,m=20.對于數(shù)值波浪水槽的入射波采用二階Stokes 波，勢函數(shù)和自由表面高程函數(shù)定義為

式中，波數(shù)k=π，波長λ =2，圓頻率ω =5．54，波高H=0．04.

由于滿足Laplace 方程的非奇異半解析基函數(shù)為平移不變調(diào)和函數(shù)（式（10）），式中的形狀參數(shù)對數(shù)值準確性至關(guān)重要，為了選取合適的形狀參數(shù)c，使用以下公式確定自由表面高程的均方根誤差：

式中，Nnts是自由液面xj∈[λ,2λ]的節(jié)點個數(shù)，ηA和 ηN分別是自由表面高程解析解（式（22））和數(shù)值結(jié)果.

此外，為了求解自由液面上速度勢 φ及其相應(yīng)偏導數(shù)值? φ/?x和?2φ/?x2的均方根誤差，使用以下公式確定均方根誤差：

式中，δA和 δN分別是自由液面φ，? φ/?x或者?2φ/?x2的解析解和數(shù)值結(jié)果.

首先，形狀參數(shù)c是影響局部邊界節(jié)點法求解精度的首要因素，因此確定最佳形狀參數(shù)是至關(guān)重要的.表1給出了在不同形狀參數(shù)下的均方根誤差.從表中可以看出，c∈[0．1,3]的均方根誤差都是在誤差允許范圍之內(nèi)的，并且c=1時 的均方根誤差最小，故本文取c=1作為最佳形狀參數(shù).

表1 不同形狀參數(shù)下二階Stokes 波的均方根誤差Table 1 The RMSEs for 2nd-order Stokes waves with different shape parameters

其次，節(jié)點總數(shù)和鄰近點數(shù)是影響計算精度、稀疏系統(tǒng)內(nèi)存占用和計算時間的主要因素.因此，表2和表3給出了不同總點數(shù)和鄰近點數(shù)下，局部邊界節(jié)點法與廣義有限差分法分別求解數(shù)值波浪水槽的均方根誤差.從表中可以看出，當總點數(shù)和鄰近點數(shù)較多時，局部邊界節(jié)點法與廣義有限差分法的誤差在同一數(shù)量級，但局部邊界節(jié)點法在具有較少總點數(shù)或者較少鄰近點數(shù)時也是可以達到較低誤差的.如表3所示，對于局部邊界節(jié)點法而言，隨著鄰近點數(shù)的增加，誤差先減小后增大.并且當鄰近點數(shù)m=20時，均方根誤差最小.所以求解數(shù)值波浪水槽時，取鄰近點數(shù)為20 左右均能保持較高的精度水平.綜上可知，相比于廣義有限差分法，局部邊界節(jié)點法可以通過取更少的總點數(shù)和鄰近點數(shù)，使得其稀疏矩陣需要更少的內(nèi)存和計算時間.

表2 不同總點數(shù)下LBKM 與GFDM 的均方根誤差Table 2 The RMSE1s of LBKM and GFDM under different total numbers of nodes

表3 不同鄰近點數(shù)下LBKM 與GFDM 的均方根誤差Table 3 The RMSE1s of LBKM and GFDM under different numbers of nearest nodes

此外，表4給出了不同總點數(shù)下自由液面上速度勢 φ及其相應(yīng)偏導數(shù)值? φ/?x和?2φ/?x2的均方根誤差.由表4可知，速度勢φ 及其相應(yīng)偏導數(shù)值? φ/?x和?2φ/?x2的均方根誤差均隨總點數(shù)增加而減小.同時，類似于傳統(tǒng)配點技術(shù)，局部邊界節(jié)點法計算得到的數(shù)值結(jié)果隨著偏導數(shù)階數(shù)的升高，誤差逐漸增大，不過其對應(yīng)的二階偏導數(shù)值誤差均低于4 ×10?2.

表4 不同總點數(shù)下φ，? φ/?x 和? 2φ/?x2的均方根誤差Table 4 The RMSE2s of φ , ? φ/?x and ? 2φ/?x2 under different total numbers of nodes

從上述討論可以看出，局部邊界節(jié)點法能夠準確、穩(wěn)定地求解數(shù)值水槽問題.圖2展示了不同總點數(shù)、不同時間間隔和不同鄰近點數(shù)下所求得的自由液面在x=4處的高程演化曲線與二階Stokes 波解析解的比較.從圖2(a)～(c)的比較結(jié)果中可以看出，除了最初幾個周期外，局部邊界節(jié)點法數(shù)值解與二階Stokes 波解析解基本吻合.由于上游入射邊界條件引入了斜坡函數(shù)，所以局部邊界節(jié)點法數(shù)值解和解析解在最初存在一定差異是合理的.圖2(a)～(c)的比較結(jié)果驗證了所提局部邊界節(jié)點法解決方案的準確性、穩(wěn)定性和一致性.

圖2 在x =4處自由液面高程演化圖：(a)不同總節(jié)點數(shù)；(b)不同時間間隔；(c)不同鄰近點數(shù)Fig.2 Evolution of free-surface elevation at x =4: (a)different total numbers of nodes; (b)different time increments;(c)different numbers of nearest nodes

為了避免波浪在下游邊界的反射現(xiàn)象，我們在下游邊界設(shè)置了海綿層達到吸收反射波的效果，并且海綿層的厚度設(shè)置為2.為了驗證海綿層的有效性，圖3展示了自由液面在4 個不同時刻的輪廓圖，從圖中可以很明顯地看出波浪呈周期性變化，但當波浪到達海綿層以后，波高呈現(xiàn)逐漸減小的趨勢，直至為零，所以波浪在下游邊界的反射現(xiàn)象基本上不可見，驗證了下游海綿層的有效性.

圖3 自由表面在4 個不同時刻的輪廓Fig.3 Profiles of free surface along the flume at 4 specific moments

2.2 波浪通過梯形防波堤的模擬

含有梯形防波堤的數(shù)值波浪水槽寬度b=30、 水深h=0．4.水下梯形防波堤的尺寸和位置如圖4所示.在此算例中，x=0處 的入射條件采用二階Stokes 波，勢函數(shù)如式（21）所示，其中波長 λ 、波周期T和波高H分別為3.693，2 和0.02.總節(jié)點數(shù)N=20 415，時間間隔 Δt=0．005，鄰近點數(shù)m=20，總時間為12T.在本例中，取Δx=0．025將節(jié)點沿x方向均勻分布.對于每個垂直方向，從實心底部到靜水表面均勻放置15 個節(jié)點，梯形防波堤的節(jié)點布置如圖5所示.

圖4 梯形防波堤數(shù)值波浪水槽的示意圖Fig.4 The schematic diagram of the numerical wave flume with a trapezoidal submerged obstacle

圖5 梯形防波堤數(shù)值波浪水槽的布點圖Fig.5 Distributions of nodes of the numerical wave flume with a trapezoidal submerged obstacle

取自由液面上的5 個點用來記錄在最后兩個模擬周期(10T～12T)的高程演化曲線，其中5 個記錄點的具體位置分別為x=5．7,x=10．5,x=12．5,x=15．7,x=17．3.5 個記錄點的高程演化曲線如圖6(a)～(e)所示，并與Ohyama 等的實驗結(jié)果[19]、廣義有限差分法結(jié)果（鄰近點數(shù)m=40）進行比較，可以觀察到局部邊界節(jié)點法結(jié)果和已有結(jié)果具有良好的一致性.從圖中可以看出，在x=5．7處，由于波浪還未與梯形防波堤發(fā)生相互作用，波峰波谷仍然呈現(xiàn)周期性變化；當波浪與梯形防波堤發(fā)生作用后，波峰變得比較陡峭，而波谷變得相對比較平坦，在x=12．5處 表現(xiàn)最為明顯，這是因為在較淺的流體區(qū)域內(nèi)，波的非線性進一步增加.此外，自由液面從x=6到x=17在 4 個特定時刻（t=10T,t=10．25T,t=10．5T,t=10．75T）的表面輪廓如圖7(a)～(d)所示.由圖7可以觀察到非線性水波與水下障礙物之間的相互作用.通過比較結(jié)果證明，局部邊界節(jié)點法能夠在較小離散區(qū)域內(nèi)準確模擬含有水下障礙物的波浪傳播問題.

圖6 最后兩個周期的高程演化Fig.6 Elevation evolution of the last 2 periods

圖7 自由液面從 x =6 到 x =17在特定時刻的表面輪廓Fig.7 Free surface profiles from x =6 to x =17 at specific moments

3 主要結(jié)論

本文首次采用局部邊界節(jié)點法結(jié)合二階Runge-Kutta 法求解數(shù)值波浪水槽.在數(shù)值波浪水槽中，首先通過2.1 小節(jié)的算例確定了形狀參數(shù)的取值范圍為c∈[0．1,3]，得到合理形狀參數(shù)c=1；并且說明局部邊界節(jié)點法能夠在較少的總節(jié)點數(shù)和較小的局部子域內(nèi)精確求解數(shù)值波浪水槽，更加適合大規(guī)模流體力學問題的計算，為大規(guī)模流體問題的求解提供了新思路；同時驗證了局部邊界節(jié)點法求解的準確性、穩(wěn)定性和一致性.其次通過2.2 小節(jié)的算例模擬了波浪與近海岸防波堤的相互作用.結(jié)果表明，當波浪與梯形防波堤發(fā)生作用后，波峰變得比較陡峭，而波谷變得相對比較平坦，這是由于在較淺的流體區(qū)域內(nèi)波的非線性進一步增加的結(jié)果.通過比較可以看出，使用局部邊界節(jié)點法能夠在較小離散區(qū)域內(nèi)精確地模擬波浪與近海岸防波堤的相互作用，從而為近海岸防波堤的相關(guān)研究和設(shè)計提供了數(shù)值參考.