水槽動(dòng)力特性數(shù)值模擬的新型局部無(wú)網(wǎng)格配點(diǎn)法*

曾維鴻，傅卓佳，湯卓超

（河海大學(xué) 力學(xué)與材料學(xué)院，南京 211100）

引 言

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

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

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

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

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

1 數(shù)學(xué)模型

1.1 控制方程

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

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

1.2 邊界條件

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

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

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

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

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

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

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

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

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

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

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

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

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

式（9）可以改寫(xiě)成如下矩陣形式：

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

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

式中

因此，式（13）可以改寫(xiě)為

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

同時(shí)，點(diǎn)x(i)處的一階和二階導(dǎo)數(shù)可以表示為

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

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

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

式中，Δt為時(shí)間增量.當(dāng)離散式（6）時(shí)，f與式（6）右端項(xiàng)相同，其中泛函V表示速度勢(shì)φ；類(lèi)似地，當(dāng)離散式（7）時(shí)，f與式（7）右端項(xiàng)相同，其中泛函V表示自由表面高程η.

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

2.1 二階Stokes 波傳播

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 在x =4處自由液面高程演化圖：(a)不同總節(jié)點(diǎn)數(shù)；(b)不同時(shí)間間隔；(c)不同鄰近點(diǎn)數(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)象，我們?cè)谙掠芜吔缭O(shè)置了海綿層達(dá)到吸收反射波的效果，并且海綿層的厚度設(shè)置為2.為了驗(yàn)證海綿層的有效性，圖3展示了自由液面在4 個(gè)不同時(shí)刻的輪廓圖，從圖中可以很明顯地看出波浪呈周期性變化，但當(dāng)波浪到達(dá)海綿層以后，波高呈現(xiàn)逐漸減小的趨勢(shì)，直至為零，所以波浪在下游邊界的反射現(xiàn)象基本上不可見(jiàn)，驗(yàn)證了下游海綿層的有效性.

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

2.2 波浪通過(guò)梯形防波堤的模擬

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

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

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

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

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

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

3 主要結(jié)論

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