去奇異化Rankine 源分布面板法自由面數(shù)值模擬

韓永浩, 羅志強(qiáng)

(昆明理工大學(xué)理學(xué)院，昆明 650500)

1 引言

流體力學(xué)中的數(shù)值模擬被廣泛地應(yīng)用于海洋運(yùn)輸、航天科技和核工程等領(lǐng)域中，二維近水面波與水翼之間的相互作用的數(shù)值模擬是近些年來(lái)流體水動(dòng)力學(xué)研究的重要課題.通過(guò)對(duì)二維水下振動(dòng)體的數(shù)值模擬，探索近水面波與水翼之間的相互作用下的振動(dòng)體周圍壓力以及自由面波等復(fù)雜流場(chǎng)內(nèi)部機(jī)制，對(duì)近水面探索器以及水下航行器等的研究具有重大意義.

在過(guò)去的幾十年中，有大量的數(shù)值模擬、理論分析和實(shí)驗(yàn)研究的方法研究近水面波與水翼之間的相互作用的問(wèn)題.Wehausen[1]、Frank[2]利用格林函數(shù)數(shù)值求解自由面波與水翼相互作用的壓力和波高.Lamb[3]建立源面板法，求解自由面非線性波.Vugts[4]計(jì)算出自由面搖擺，起伏和滾動(dòng)圓柱三種情況下的水動(dòng)力系數(shù).Cao 等[5]運(yùn)用去奇異化的邊界方法求解出三維勢(shì)流問(wèn)題.Beck[6]運(yùn)用時(shí)域面板法求解出了浮體在無(wú)限深度下波體相互作用下的水動(dòng)力系數(shù).自由面邊界條件對(duì)于處理近水域振動(dòng)體和波之間相互作用的問(wèn)題時(shí)非常重要，如何進(jìn)行邊界條件的處理是當(dāng)時(shí)一大難題.Dawson[7]在1977 年提出了Ranine 源方法，使用Rankine 源法滿足輻射和自由面邊界條件，用于描述勻速前進(jìn)的物體的穩(wěn)定的流動(dòng)，并獲得了成功.Beck 等[8]運(yùn)用Rankine 源法，把自由面邊界積分和圓柱體邊界積分兩個(gè)方程分開(kāi)表達(dá)求解，數(shù)值模擬了完全非線性波.由于在自由面的邊界積分涉及到無(wú)限自由面，Lee[9]將自由面區(qū)域截?cái)酁橛邢迏^(qū)域進(jìn)行積分求解，解決了表面穿透體的非線性輻射問(wèn)題.Lee[10]考慮孤立源求解完全非線性流體動(dòng)力問(wèn)題.Zhang 和Beck[11,12]考慮在有限的自由面區(qū)域均勻放置Rankine 源的方法，求解出了二維漂浮物大振幅下的問(wèn)題.Feng 等[13-16]考慮在自由面Rankine 源稀疏放置，數(shù)值模擬了水面圓柱作垂直方向運(yùn)動(dòng)的水動(dòng)力學(xué)性態(tài).

綜上所述，以上文獻(xiàn)考慮各種數(shù)值方法數(shù)值求解自由面波高與水翼表面壓力的水動(dòng)力學(xué)性態(tài).很少文獻(xiàn)考慮水下物體作垂直運(yùn)動(dòng)時(shí)波體之間的水動(dòng)力學(xué)性態(tài).由于自由面波的運(yùn)動(dòng)受到水下運(yùn)動(dòng)物體和水下源強(qiáng)度的雙重影響，很難構(gòu)建直接關(guān)系的表達(dá)式，因而給自由面波的求解造成了困難.鑒于此，本文構(gòu)建了水下圓柱作垂直運(yùn)動(dòng)下，在圓柱體和自由面同時(shí)放置Rankine 源，巧妙放置Rankine 源的具體位置，使得源的分布符合流體運(yùn)動(dòng)的規(guī)律，同時(shí)巧妙地使用了自由面插值的方法構(gòu)建自由面波波高與圓柱運(yùn)動(dòng)的關(guān)系式，求解了源分布下波體運(yùn)動(dòng)相互作用的水動(dòng)力學(xué)性態(tài).

本文第2 部分為模型的建立與邊界積分方程的推導(dǎo)，第3 部分為邊界積分方程的離散與處理，第4 部分為數(shù)值實(shí)驗(yàn)?zāi)M，最后給出了本文的結(jié)論.

2 數(shù)學(xué)模型以及積分方程的推導(dǎo)

2.1 物理模型

假設(shè)在不可壓、無(wú)粘和無(wú)旋的理想流體中，考慮二維無(wú)限水深的水面下振動(dòng)圓柱體水翼，如圖圖1 所示.當(dāng)Rankine 源Q 點(diǎn)分布在流體區(qū)域V 中，設(shè)SQ是以Q 為圓心，ε 為半徑的圓的邊界.流體區(qū)域V 由自由面邊界Sf、圓柱體邊界Sb、源點(diǎn)Q 所圍成的區(qū)域邊界SQ和無(wú)窮遠(yuǎn)處海底邊界圍成，其中Sc為平靜水面邊界，R 為振動(dòng)圓柱體半徑，h 為振動(dòng)體浸沒(méi)深度.

圖1 物理模型

圓柱體垂直運(yùn)動(dòng)速度為

U =(0,a sin ωt),

其中U 為圓柱體垂直運(yùn)動(dòng)速度，a 為振動(dòng)的振幅，ω 為振動(dòng)頻率.

流體區(qū)域內(nèi)勢(shì)流φ 滿足?2φ=0.圓柱體與液體的固液邊界條件為

其中nB為圓柱體表面外法向量.

自由面運(yùn)動(dòng)學(xué)邊界條件為

其中y =η(x,t)為自由面波高.

自由面動(dòng)力學(xué)邊界條件可以表示為

其中g(shù) 為重力加速度.

2.2 源分布邊界積分方程

設(shè)φ 和φ1是流體區(qū)域V 的兩種速度勢(shì)函數(shù)，分別滿足Laplace 方程

?2φ=0, ?2φ1=0.

φ 和φ1在流體區(qū)域V 中是有限單值函數(shù).對(duì)應(yīng)的速度勢(shì)關(guān)于x 和y 的導(dǎo)數(shù)

設(shè)S 是流體區(qū)域V 的所有邊界之和，其中l(wèi), w 是邊界S 上外法向量n 在x, y 正方向上的單位分量，得

由格林公式

得

取φ1=ln r =ln|P -Q|，其中Q 點(diǎn)是源點(diǎn)，P 點(diǎn)為流體區(qū)域V 內(nèi)部的點(diǎn).

當(dāng)源點(diǎn)Q 在流體域V 內(nèi)部，如圖2(a)，P 取Q 點(diǎn)處時(shí)ln r 無(wú)限大，積分區(qū)域?yàn)閂 -VQ.流體區(qū)域V 的邊界為

(5)式改寫(xiě)為

圖2 源點(diǎn)Q 位置以及內(nèi)外流場(chǎng)剖面

根據(jù)

(6)式化為

當(dāng)ε →0 時(shí)

當(dāng)源點(diǎn)Q 在流體區(qū)域的外部(圖2(a)中Q1所處位置)，流體區(qū)域內(nèi)的點(diǎn)P 到源點(diǎn)Q 的函數(shù)ln r 取值有限.流體區(qū)域V 的邊界為

(5)式改寫(xiě)為

設(shè)源點(diǎn)Q 影響的區(qū)域?yàn)榈谝粎^(qū)域，源點(diǎn)Q 在第一區(qū)域的內(nèi)部，φ 是第一區(qū)域速度勢(shì)，φ1是流體區(qū)域V 剩余區(qū)域的速度勢(shì)，如圖2(b)，由(7),(8)式推出對(duì)應(yīng)的邊界積分

由(9)和(10)式可以推出

令

因?yàn)棣?φ1，且速度勢(shì)φ, φ1關(guān)于邊界S1的法向?qū)?shù)不連續(xù)，得到圓柱邊界積分方程

3 Rankine 源分布邊界積分方程離散

3.1 邊界積分方程

3.2 積分方程的離散化

圖3 t=0 時(shí)刻自由面和圓柱面面板剖分圖

3.3 時(shí)間域中源點(diǎn)、控制點(diǎn)、切向量

在t = 0 自由面上放置源點(diǎn)時(shí)，由于流體區(qū)域是關(guān)于振動(dòng)體對(duì)稱的，因此，我們只討論y 軸右邊的自由面源點(diǎn)的放置方法.在本文中選取Nf+1=181 個(gè)源點(diǎn)，所以在右邊自由面有Nf/2=90 個(gè)源點(diǎn)，70 個(gè)源點(diǎn)在內(nèi)域，20 個(gè)源點(diǎn)在外域.內(nèi)域源點(diǎn)放置方法

外域源點(diǎn)放置方法

其中Lb為圓柱體面板長(zhǎng)度，αj設(shè)置如表1.

表1 αj 的參數(shù)

由于振動(dòng)體圓柱完全浸沒(méi)在自由面水下，振動(dòng)體面板進(jìn)行均勻分割，因此振動(dòng)體圓柱體的面板長(zhǎng)度并不隨時(shí)間的變化進(jìn)行變化，同理物體表面單個(gè)面板的切向量亦不隨時(shí)間作變化.設(shè)圓柱體表面Sb的源點(diǎn)，面板長(zhǎng)度，面板切向量分別為

隨著振動(dòng)時(shí)間的變化，自由面隨圓柱體振動(dòng)進(jìn)行變化.設(shè)自由面Sf的源點(diǎn)，面板長(zhǎng)度，面板切向量分別為

自由面邊界積分方程

同時(shí)，圓柱體邊界滿足積分方程

通過(guò)邊界條件(1)和(15)式邊界條件結(jié)合可以得到圓柱體邊界積分方程

(14)式和(16)式化為下面系列代數(shù)方程

其中σi為

bi為

bi=Uni, i=1,2,··· ,Nb, bi=φi, i=1+Nb,··· ,N,

其中影響系數(shù)ai,j為

壓力P 的表達(dá)式為

自由面波高η 表達(dá)式為

4 數(shù)值模擬

自由面波關(guān)于y 軸對(duì)稱，下列數(shù)值模擬中，取振動(dòng)圓柱體半徑R = 10 m，圓柱體面板數(shù)Nb= 40，自由面面板數(shù)Nf= 180，時(shí)間步長(zhǎng)△t = 0.001 s，浸沒(méi)深度h，圓柱振動(dòng)頻率ω，圓柱振動(dòng)振幅a，振動(dòng)時(shí)間t，自由面波畫(huà)圖時(shí)選取原點(diǎn)處源點(diǎn).

時(shí)域函數(shù)迭代的初始自由面邊界條件

4.1 數(shù)值驗(yàn)證

為了驗(yàn)證文中采用的數(shù)值方法的準(zhǔn)確性和有效性，本文采用自由面波高與其他文獻(xiàn)中的數(shù)值解進(jìn)行比較.

圖4 采用圓柱體振動(dòng)頻率ω = 1.0，振幅a = 1m.圖中表示自由面波的數(shù)值解與文獻(xiàn)[16]的數(shù)值解比較圖，從圖中可以發(fā)現(xiàn)結(jié)果吻合得很好，表明本文采用的數(shù)值方法是準(zhǔn)確有效的.

圖4 自由面波高

4.2 自由面波和圓柱體壓力數(shù)值模擬

首先，圖5(a)研究了當(dāng)浸沒(méi)深度h/R = 1.5 時(shí)，不同的圓柱體振動(dòng)頻率下自由面波高變化規(guī)律，從圖5(a)得出當(dāng)振動(dòng)頻率增大的時(shí)候，自由面波高也相應(yīng)的增大.圖5(b)中取ω = 0.5，研究不同的浸沒(méi)深度，得出浸沒(méi)深度越大，自由面波高越小.對(duì)比圖5(b)和圖5(c)在不同浸沒(méi)深度下的自由面波可以看到，振動(dòng)頻率越大，自由面波高也大，說(shuō)明這種論證不僅局限于一定的浸沒(méi)深度，對(duì)不同的浸沒(méi)深度也有這種效果.對(duì)比圖5(c)和圖5(d)，更加論證了上述浸沒(méi)深度越小，自由面波高越大的觀點(diǎn).圖5(f)取浸沒(méi)深度h/R =1.5 一定時(shí)，數(shù)值模擬了不同振動(dòng)振幅下，振動(dòng)體對(duì)自由面波的影響情況，從圖中可以看出，振動(dòng)體振幅越大，自由面波波高越大，影響效果越明顯.

圖5 不同振動(dòng)參數(shù)下自由面波高變化規(guī)律

最后數(shù)值模擬了不同參數(shù)下，振動(dòng)體圓柱體表面壓力情況.在圖6(a)，取浸沒(méi)深度h/R = 1.2，數(shù)值模擬了不同時(shí)刻下的圓柱體表面壓力情況，可以得到不同時(shí)刻下，圓柱體的壓力隨圓柱體的振動(dòng)在進(jìn)行變化，而不是固定不變的.圖6(b)中，取浸沒(méi)深度h/R=1.5，對(duì)比圖6(a)和圖6(b)可以得到，浸沒(méi)深度的越小，振動(dòng)圓柱體表面壓力越大.為了更形象的說(shuō)明上述觀點(diǎn)，在圖6(c)中，取固定振動(dòng)時(shí)刻t = 10s，在圖6(c)中，可以更清晰的論證上述觀點(diǎn).圖6(d)模擬了不同振動(dòng)振幅下，圓柱體周圍壓力變化情況，從圖6(d)中可以得出.振動(dòng)振幅越大，振動(dòng)圓柱體周圍壓力越大.

圖6 不同振動(dòng)參數(shù)下圓柱體表面壓力變化規(guī)律

5 結(jié)論

本文基于源面板法，在二維無(wú)旋、無(wú)粘、不可壓縮的流體中，數(shù)值模擬了水下振動(dòng)體振動(dòng)對(duì)自由面波高和圓柱體壓力的影響.通過(guò)本文的數(shù)值模擬結(jié)果可知，圓柱體振動(dòng)自由面波高和壓力的變化與圓柱體浸沒(méi)深度、圓柱體振動(dòng)頻率以及振動(dòng)體振動(dòng)振幅有關(guān).隨振動(dòng)體振動(dòng)頻率或振動(dòng)振幅的增大，自由面波高也隨之增大，減小振動(dòng)體浸沒(méi)深度，自由面波高也逐漸增加.同時(shí)還得到，振動(dòng)體振動(dòng)振幅的增大也導(dǎo)致圓柱體表面壓力的增大，隨振動(dòng)體浸沒(méi)深度的增大，振動(dòng)體周圍壓力反而減小.