， ， ， ， ，

(1. 大連理工大學(xué) 船舶工程學(xué)院， 遼寧 大連， 116024； 2. 大連海洋大學(xué) 海洋與土木工程學(xué)院， 遼寧 大連， 116023)

0 引 言

海洋結(jié)構(gòu)物是海洋資源開(kāi)發(fā)長(zhǎng)期作業(yè)的依托平臺(tái)，開(kāi)發(fā)過(guò)程中往往伴隨著惡劣的海洋環(huán)境條件，引起結(jié)構(gòu)物的不利荷載和劇烈運(yùn)動(dòng)[1]，在極端情況下，甚至?xí)鸾Y(jié)構(gòu)物整體傾覆，導(dǎo)致重大的工程事故。但是，傳統(tǒng)水動(dòng)力軟件如AQWA[2]、SESAM等都采用線性波浪假定，忽略了波浪的非線性影響，無(wú)法滿足極限波浪情況下海洋結(jié)構(gòu)的設(shè)計(jì)需要。因此，有必要對(duì)非線性波浪與結(jié)構(gòu)物的作用問(wèn)題做進(jìn)一步研究，在設(shè)計(jì)階段準(zhǔn)確地預(yù)測(cè)出海洋結(jié)構(gòu)物在服役期間可能受到海洋環(huán)境的最大沖擊，避免海洋工程事故的發(fā)生。

目前，關(guān)于完全非線性波浪與結(jié)構(gòu)物相互作用的研究有2種不同的研究方向：一種方法統(tǒng)稱數(shù)值波浪水槽技術(shù)[3]，模擬物理試驗(yàn)水槽技術(shù)建立數(shù)值波浪水槽，結(jié)構(gòu)物放置于數(shù)值水槽中央，兩端使用側(cè)壁，前后分別設(shè)置造波與消波邊界，模擬整個(gè)造波、波物作用與消波的過(guò)程。另一種方法則基于入散射勢(shì)分離技術(shù)開(kāi)發(fā)的開(kāi)敞水域技術(shù)[4]，將速度勢(shì)分離成入射勢(shì)和散射勢(shì)，對(duì)散射勢(shì)建立積分方程，在一個(gè)開(kāi)敞圓域內(nèi)建立波浪與結(jié)構(gòu)物相互作用的模型。由于開(kāi)敞水域技術(shù)更符合工程實(shí)際情況，且具有便于造波、消波以及計(jì)算量少、計(jì)算效率高等優(yōu)點(diǎn)，逐漸被學(xué)者研究和采用。2013年，ZHOU等[4]采用開(kāi)敞水域完全非線性技術(shù)對(duì)漂浮圓柱體的繞射、受迫振動(dòng)以及Ringing現(xiàn)象進(jìn)行深入的研究，揭示非線性作用的機(jī)理。2014年，BAI等[5]對(duì)完全浸沒(méi)水中的自由及系泊的柱狀結(jié)構(gòu)物進(jìn)行完全非線性模擬。2015年，F(xiàn)ENG等[6]利用完全非線性數(shù)值模型對(duì)并列雙箱在波浪作用下的水動(dòng)力共振問(wèn)題進(jìn)行研究。2016年，HANNAN等[7]利用完全非線性水波理論對(duì)完全浸沒(méi)水中的起重船的有效工作載荷進(jìn)行數(shù)值模擬分析。但是，上述研究?jī)H限于考慮波浪與圓柱以及方箱等簡(jiǎn)單幾何形狀的物體作用。這主要是由于完全非線性波物作用理論要求自由水面與物面條件均在瞬時(shí)滿足，計(jì)算中需要在每一時(shí)刻對(duì)自由水面以及物面進(jìn)行更新。當(dāng)考慮波浪對(duì)復(fù)雜曲面結(jié)構(gòu)的作用時(shí)，復(fù)雜幾何形狀的數(shù)值計(jì)算極易導(dǎo)致數(shù)值不穩(wěn)定的情況發(fā)生。因此，建立精確自由水面及物面網(wǎng)格的更新方法，準(zhǔn)確模擬復(fù)雜曲面波物運(yùn)動(dòng)過(guò)程，是完全非線性波物作用模型應(yīng)用于工業(yè)設(shè)計(jì)所亟需解決的關(guān)鍵性問(wèn)題。

本文采用開(kāi)敞水域數(shù)值模擬技術(shù)以及高階邊界元法[8]，在開(kāi)敞圓域內(nèi)建立一個(gè)完全非線性波浪與海洋結(jié)構(gòu)物作用的三維數(shù)值分析模型。通過(guò)完全非線性波浪對(duì)直立圓柱作用的高階波浪力進(jìn)行模擬分析，驗(yàn)證本文模型的正確性。為驗(yàn)證本文模型對(duì)復(fù)雜曲面問(wèn)題的適用性，對(duì)非線性波浪與Wigley船的相互作用問(wèn)題進(jìn)行模擬，并通過(guò)與JOURNEE[9]的物理模型試驗(yàn)對(duì)比，驗(yàn)證模擬方法的可靠性。

1 控制方程與邊界條件

圖1 計(jì)算模型示例

波浪與結(jié)構(gòu)物作用問(wèn)題的計(jì)算模型如圖1所示(圖中：SD為水底固面；SB為物面邊界；SF為自由水面)。定義2組坐標(biāo)系O-xyz和O-x′y′z′分別描述整個(gè)流體域和物體的運(yùn)動(dòng)，法向量定義為出流體方向?yàn)檎?/p>

在勢(shì)流理論的假定下，以三維拉普拉斯方程[10]為基本控制方程，在時(shí)域內(nèi)建立完全非線性波浪與浮體間作用的數(shù)學(xué)模型，并且瞬時(shí)自由水面條件和物面條件均是非線性的。根據(jù)入散射勢(shì)分離技術(shù)，將總速度勢(shì)和波面(φ，η)分解成入射勢(shì)(φI，ηI)和散射勢(shì)(φS，ηS)，代入控制方程和邊界條件并整理，得

2φS=0 (在計(jì)算域D中) (1)

式中：U為結(jié)構(gòu)物的平動(dòng)速度；Ω為結(jié)構(gòu)物的轉(zhuǎn)動(dòng)速度；t為時(shí)間；n為法向量；rb為任一點(diǎn)到轉(zhuǎn)動(dòng)中心的距離。

采用五階Stokes波理論計(jì)算入射勢(shì)的有關(guān)項(xiàng)，并乘以緩沖函數(shù)來(lái)避免初始效應(yīng)。為保證輻射波浪向外傳播而不反射回來(lái)，在計(jì)算域外圍布置一層人工阻尼層。

2 模型建立

2.1 邊界積分方程建立

采用Rankine源格林函數(shù)[11]，在整個(gè)計(jì)算域內(nèi)應(yīng)用格林第二定理，將散射勢(shì)φS的邊值問(wèn)題轉(zhuǎn)化為邊界積分方程問(wèn)題：

(7)

式中：c為固角系數(shù)；G為格林函數(shù)。

通過(guò)高階邊界元法將積分方程進(jìn)行離散，可得

(ξ,ζ)|(8)

式中：Hk為權(quán)系數(shù)；N為離散節(jié)點(diǎn)總數(shù)；m為任一離散單元；M為離散單元總數(shù)；K為任一單元上的節(jié)點(diǎn)總數(shù)；j為離散單元上節(jié)點(diǎn)；J為雅克比函數(shù)；ξ、ζ分別為任一節(jié)點(diǎn)在參數(shù)坐標(biāo)系下的橫、縱坐標(biāo)值。

其中，Rankine源格林函數(shù)為

2.2 復(fù)雜船體曲面的網(wǎng)格更新

如何準(zhǔn)確地追蹤自由水面上流體質(zhì)點(diǎn)和復(fù)雜船體曲面上濕表面的變化情況，并根據(jù)相應(yīng)的變化對(duì)水面和物面上的網(wǎng)格節(jié)點(diǎn)進(jìn)行更新與重構(gòu)，是完全非線性波浪與三維船體相互作用模擬的重點(diǎn)也是難點(diǎn)。

本文采用混合歐拉-拉格朗日法(Mixed Euler-Lagrange method，MEL)追蹤自由水面流體質(zhì)點(diǎn)，采用跟隨式網(wǎng)格拓?fù)浼夹g(shù)保證網(wǎng)格單元的高質(zhì)量。自由水面節(jié)點(diǎn)位置可以用水平坐標(biāo)(x,y)和垂向坐標(biāo)z(η)來(lái)表示。在每一時(shí)間步內(nèi)，首先通過(guò)MEL對(duì)各節(jié)點(diǎn)位置坐標(biāo)進(jìn)行更新。為保持在水平投影面上各節(jié)點(diǎn)與運(yùn)動(dòng)物體的水平相對(duì)位置不變，通過(guò)初始時(shí)刻的網(wǎng)格節(jié)點(diǎn)水平坐標(biāo)和剛體水平位移得到新的網(wǎng)格節(jié)點(diǎn)水平坐標(biāo)

Δxi(i=1,2)(10)

最后，通過(guò)插值技術(shù)，根據(jù)新的節(jié)點(diǎn)水平坐標(biāo)，在最開(kāi)始運(yùn)用MEL得到的水面節(jié)點(diǎn)坐標(biāo)中，插值得到新節(jié)點(diǎn)處的波面高程和速度勢(shì)。首先找到距離目標(biāo)點(diǎn)最近的節(jié)點(diǎn)，在其所在的每個(gè)單元中有

(11)

式中：h為形函數(shù)。

因?yàn)閔只與參數(shù)坐標(biāo)(ξ,η)有關(guān)，將上式改寫為

(12)

采用平面二分法插值可以求得參數(shù)坐標(biāo)(ξ,η)，進(jìn)而通過(guò)型函數(shù)得到新節(jié)點(diǎn)的波面高程和速度勢(shì)，得到新的相對(duì)均勻的自由水面單元節(jié)點(diǎn)參數(shù)。處于自由水面和物面相交的水線上的節(jié)點(diǎn)具有特殊性(同屬于水面和物面)。為了保證經(jīng)過(guò)自由水面條件更新后的水線節(jié)點(diǎn)仍處于物面上，需要對(duì)其進(jìn)行特殊處理。

圖2 水線節(jié)點(diǎn)處理示例

圖3 船體表面網(wǎng)格重構(gòu)示例

圖2為空間坐標(biāo)系和隨體坐標(biāo)系下的船體橫剖面圖。點(diǎn)P0(x0,y0,z0)位于水線上，經(jīng)過(guò)自由水面條件更新后，與船體表面分離，需要將其移至船體表面上的目標(biāo)點(diǎn)P(x,y,z)處。首先，通過(guò)坐標(biāo)轉(zhuǎn)換矩陣得到點(diǎn)P0在隨體坐標(biāo)系中對(duì)應(yīng)的坐標(biāo)點(diǎn)P0′(x0′,y0′,z0′)；然后，根據(jù)船體橫剖曲線關(guān)系(Wigley船體函數(shù)關(guān)系或者船體橫剖線圖)將P0′強(qiáng)行移動(dòng)到船體表面上，得到點(diǎn)P1′(x0′,y1′,z0′)；最后，運(yùn)用轉(zhuǎn)換矩陣，得到空間坐標(biāo)系中船體表面上的點(diǎn)P1(x1,y1,z1)。由圖2可以看出：點(diǎn)P1與目標(biāo)點(diǎn)P的位置非常接近。因?yàn)槟M中船體在每一時(shí)間步的位移值極小，可以近似認(rèn)為點(diǎn)P1即為交線上的點(diǎn)。

船體濕表面形狀不斷變化，物面網(wǎng)格也需要重構(gòu)。船體曲面形狀的準(zhǔn)確捕捉和還原對(duì)船體在非線性波浪作用下的受力和運(yùn)動(dòng)響應(yīng)的模擬準(zhǔn)確性非常重要。在隨體動(dòng)坐標(biāo)系下，假設(shè)物面網(wǎng)格節(jié)點(diǎn)的橫坐標(biāo)x在整個(gè)模擬過(guò)程中保持不變，改變各節(jié)點(diǎn)的y軸和z軸坐標(biāo)值。如圖3所示，首先根據(jù)水線節(jié)點(diǎn)的改變量，通過(guò)彈簧近似法[4]對(duì)物面各節(jié)點(diǎn)的z坐標(biāo)進(jìn)行處理：

Δzn+1=∑kijΔzn/∑kij(13)

式中：Δz為各節(jié)點(diǎn)z坐標(biāo)的改變量；kij為彈簧剛度系數(shù)；n為迭代次數(shù)。

最后，由節(jié)點(diǎn)的x坐標(biāo)以及移動(dòng)后的z坐標(biāo)，根據(jù)船體表面的曲面函數(shù)或者通過(guò)橫剖線圖進(jìn)行插值，即可得到移動(dòng)后節(jié)點(diǎn)的y坐標(biāo)值。上述物面網(wǎng)格重構(gòu)方法對(duì)實(shí)際曲面具有較高的還原性和通用性，適用于各類復(fù)雜曲面海洋結(jié)構(gòu)物。

綜上所述，經(jīng)過(guò)自由水面網(wǎng)格的更新、水線節(jié)點(diǎn)的處理以及船體表面網(wǎng)格的重構(gòu)等3個(gè)處理過(guò)程，就可以實(shí)現(xiàn)捕捉瞬時(shí)液面和瞬時(shí)復(fù)雜船體濕表面的目標(biāo)，進(jìn)而實(shí)現(xiàn)對(duì)整個(gè)完全非線性過(guò)程的時(shí)域模擬。整個(gè)處理過(guò)程具有較強(qiáng)的適用性，理論上可以處理實(shí)際海洋工程中大多數(shù)具有復(fù)雜曲面形狀的結(jié)構(gòu)物。

2.3 ?φ/?n奇異性處理

當(dāng)計(jì)算域邊界表面(如非線性自由水面)變化較大時(shí)，格林函數(shù)的導(dǎo)數(shù)?φ/?n會(huì)產(chǎn)生柯西主值積分奇異性問(wèn)題，對(duì)完全非線性的時(shí)域模擬造成困難。?φ/?n奇異性的處理方法有2種：一種是BAI等[5]采用的間接計(jì)算法，方法簡(jiǎn)單但是需要增加額外的計(jì)算控制面，增加了數(shù)值模擬的計(jì)算量與時(shí)間；另一種是NING等[12]采用的直接計(jì)算法，可以得到一定精確度的奇異積分值，但是當(dāng)源點(diǎn)位于單元邊界中點(diǎn)時(shí)效果不是很好。

本文在第2種方法的基礎(chǔ)上進(jìn)行適當(dāng)改進(jìn)，通過(guò)單元細(xì)分、積分域分解、極坐標(biāo)變換以及級(jí)數(shù)展開(kāi)等過(guò)程將奇異積分項(xiàng)逐步轉(zhuǎn)化為常規(guī)積分，實(shí)現(xiàn)奇異項(xiàng)的消除。包含格林函數(shù)導(dǎo)數(shù)部分的積分在轉(zhuǎn)換后的局部坐標(biāo)系中可以表示為

式中：φ(x)為速度勢(shì)；Tp為參數(shù)坐標(biāo)系下的實(shí)際單元的映射；h為形函數(shù)，代表單元節(jié)點(diǎn)對(duì)單元內(nèi)任一點(diǎn)參數(shù)值的貢獻(xiàn)。

圖4 參數(shù)空間中的極坐標(biāo)系

當(dāng)源點(diǎn)與j節(jié)點(diǎn)重合時(shí)，式(14)中第j項(xiàng)的積分即為奇異積分。如圖4所示，根據(jù)與源點(diǎn)重合的節(jié)點(diǎn)位置不同，可以將單元節(jié)點(diǎn)分為角點(diǎn)和邊中點(diǎn)等2種情況。以第1種情況為例，采用直接法求解?φ/?n奇異積分。首先將參數(shù)單元分為4個(gè)小正方形單元，在小單元2、3、4中無(wú)奇異性影響，仍采用常規(guī)高斯積分。在1號(hào)小單元中，可以看作奇異性存在于α-β之間的區(qū)域內(nèi)，在這一區(qū)域內(nèi)積分可以表示為

ξdζ(15)

在剩下的區(qū)域中積分為

ξdζ(16)

式中：Tp2為單元區(qū)域減去奇異性區(qū)域后剩下的區(qū)域。

將實(shí)際區(qū)域α-β看作平面，在源點(diǎn)附近可以對(duì)空間向量進(jìn)行泰勒展開(kāi)：

代入式(16)，可以得到

(18)

式中：θ為極坐標(biāo)系下角度的參數(shù)；f(θ)為表達(dá)式的函數(shù)；R為單元邊界上點(diǎn)距源點(diǎn)的距離；α為奇異區(qū)域外緣距源點(diǎn)的距離；ρ為極坐標(biāo)系下徑向的參數(shù)。其中：

積分IS2中包含2個(gè)同樣的奇異積分可以相互抵消，采用常規(guī)高斯積分即可。

將式(17)代入式(15)，可以得到

式中：ε為奇異區(qū)域邊界半徑；A為定義的偏導(dǎo)數(shù)參量Ai(θ)=cos(θ)?xi/?ζ+sin(θ)?xi/?ζ。

式(20)等號(hào)右側(cè)第1項(xiàng)為無(wú)奇異項(xiàng)，可以通過(guò)常規(guī)高斯積分計(jì)算。式(20)等號(hào)右側(cè)第2項(xiàng)積分隨ε趨于0而趨于無(wú)窮，但當(dāng)源點(diǎn)周圍所有對(duì)應(yīng)積分相加后總和必然等于0，故可以不考慮。

通過(guò)上述單元細(xì)分、積分方程分解、極坐標(biāo)變換以及空間向量級(jí)數(shù)展開(kāi)等處理過(guò)程，可以將包含格林函數(shù)導(dǎo)數(shù)積分中的奇異積分部分消除。在保證計(jì)算精度的前提下，只需要額外增加很小的計(jì)算量，就可以得到較高的數(shù)值模擬效率，有助于對(duì)非線性波浪與復(fù)雜曲面船體相互作用的模擬模型的建立。

3 數(shù)值結(jié)果

3.1 完全非線性波浪對(duì)直立圓柱的繞射問(wèn)題

圖5 直立圓柱和自由水面初始網(wǎng)格劃分示例

本文通過(guò)與HUSEBY等[13]的試驗(yàn)結(jié)果及ZHOU等[4]的數(shù)值結(jié)果進(jìn)行對(duì)比，驗(yàn)證計(jì)算模型中波浪力計(jì)算結(jié)果的準(zhǔn)確性。采用與試驗(yàn)相同的模型進(jìn)行分析，如圖5所示，直立貫底圓柱半徑r=0.03 m，計(jì)算域水深為d=0.6 m，水域半徑取2倍波長(zhǎng)。根據(jù)對(duì)稱性，在1/2物面和自由水面上劃分網(wǎng)格，其中物面上剖分12×10個(gè)單元，在水面上剖分12×20個(gè)單元。入射勢(shì)采用5階Stokes波理論解，波數(shù)k=8.1，與試驗(yàn)值相同，僅考慮波幅A的變化對(duì)圓柱受力的影響情況。將數(shù)值試驗(yàn)結(jié)果通過(guò)傅里葉變換后與試驗(yàn)結(jié)果進(jìn)行比較。各階波浪力對(duì)比如圖6所示。

圖6 各階波浪力對(duì)比

從圖6可以看出：本文數(shù)值模型計(jì)算結(jié)果大體上與HUSEBY等[13]的試驗(yàn)結(jié)果及ZHOU等[4]的數(shù)值結(jié)果吻合，證明了數(shù)值模型的準(zhǔn)確性。在未損失太多計(jì)算精度的前提下可以模擬更大波幅波浪的入射情況。其中波浪一階力的計(jì)算結(jié)果略小于其他2種結(jié)果，但相差很小，在可接受的誤差范圍內(nèi)。隨著振幅增加波浪力變化平緩，變化趨勢(shì)與試驗(yàn)結(jié)果相同。波浪二階力作為最主要的高階力，在一定波高范圍內(nèi)，本文的計(jì)算結(jié)果與試驗(yàn)吻合良好。本文波浪三階力和四階力的結(jié)果隨著波幅A變化較大；本文五階波浪力的計(jì)算結(jié)果小于其他2種結(jié)果，但差距不大。

3.2 完全非線性波浪對(duì)Wigley船的作用問(wèn)題

為了驗(yàn)證本文數(shù)值模型對(duì)凸出水線型復(fù)雜曲面海洋結(jié)構(gòu)物在完全非線性波浪作用下的模擬能力，對(duì)零航速Wigley船在非線性波浪作用下受到的波浪力和運(yùn)動(dòng)響應(yīng)進(jìn)行模擬，并將結(jié)果與JOURNEE[9]的試驗(yàn)結(jié)果進(jìn)行對(duì)比。本文數(shù)值模擬所采用的Wigley船型的數(shù)學(xué)表達(dá)式為

圖7 Wigley船及自由水面網(wǎng)格劃分

船長(zhǎng)L=3 m，船寬B=0.3 m，吃水T=0.187 5 m，船排水量0.078 m3，重心距船底0.17 m，縱向慣性半徑Ryy=0.75。入射波采用五階Stokes波，入射角180°，波幅A=0.01 m，波長(zhǎng)λ/L取值0.5～2.5。根據(jù)對(duì)稱性，在1/2 Wigley殼體上劃分8×16個(gè)單元。計(jì)算域半徑取為R=5L，在1/2水面上劃分16×30個(gè)單元。Wigley船及自由水面網(wǎng)格劃分如圖7所示。

圖8為Wigley船受到的波浪力和力矩隨無(wú)因次波長(zhǎng)λ/L的變化情況。由圖8 a)可以看出：本文關(guān)于Wigley船的縱向波浪力的數(shù)值結(jié)果與試驗(yàn)值大部分吻合良好，只有在波長(zhǎng)較小時(shí)，兩者存在一些偏差。由圖8 b)可以看出：本文關(guān)于Wigley船垂向波浪力的模擬計(jì)算結(jié)果與試驗(yàn)值存在一定誤差，波長(zhǎng)較大時(shí)數(shù)值結(jié)果略大于試驗(yàn)值。由圖8 c)可以看出：繞y軸縱搖力矩的計(jì)算結(jié)果在波長(zhǎng)較小時(shí)與試驗(yàn)結(jié)果吻合良好，在大波長(zhǎng)時(shí)略小于試驗(yàn)值。

圖8 Wigley船波浪力和波浪力矩對(duì)比

圖9為Wigley船在波浪作用下的位移響應(yīng)與試驗(yàn)結(jié)果對(duì)比，可以看出：Wigley船的垂向位移在波長(zhǎng)較小時(shí)比試驗(yàn)值偏小且沒(méi)有反應(yīng)出試驗(yàn)值的波動(dòng)性，可能是模型試驗(yàn)精度誤差或者表面張力黏性引起的；而在波長(zhǎng)λ/L大于1.0時(shí)，數(shù)值計(jì)算結(jié)果偏大一些，但是誤差不大。本文的縱搖計(jì)算結(jié)果與試驗(yàn)結(jié)果基本吻合，除個(gè)別點(diǎn)外偏差不大，具有一定的可靠性。因此，本模型對(duì)具有非規(guī)則表面的Wigley船的時(shí)域模擬基本與試驗(yàn)值吻合良好，進(jìn)而證明了本模型的處理方法對(duì)瞬時(shí)自由水面和船體復(fù)雜濕表面的跟蹤與還原的準(zhǔn)確性，具有較高的實(shí)用價(jià)值。

圖9 Wigley船位移響應(yīng)對(duì)比

4 結(jié) 論

(1) 通過(guò)跟隨式網(wǎng)格拓?fù)浼夹g(shù)和曲面網(wǎng)格重構(gòu)技術(shù)更新網(wǎng)格節(jié)點(diǎn)坐標(biāo)，實(shí)現(xiàn)了對(duì)瞬時(shí)非線性自由水面和瞬時(shí)復(fù)雜船體濕表面的準(zhǔn)確跟蹤與還原。

(2) 通過(guò)對(duì)圓柱繞射問(wèn)題高階波浪力的模擬測(cè)試，證明了模型對(duì)瞬時(shí)自由表面變化捕捉的準(zhǔn)確性以及對(duì)流體作用力模擬的準(zhǔn)確性。

(3) 通過(guò)對(duì)零航速Wigley船完全非線性運(yùn)動(dòng)響應(yīng)的模擬，驗(yàn)證了本文數(shù)值模型對(duì)復(fù)雜船體濕表面的還原能力以及在非線性波浪作用下運(yùn)動(dòng)響應(yīng)的模擬能力。