劉夢超 ，劉延俊 ，，薛鋼 ，吳瀚崚

1山東大學(xué)海洋研究院，山東濟(jì)南250100

2山東大學(xué)機(jī)械工程學(xué)院，山東濟(jì)南250061

3山東大學(xué)高效清潔機(jī)械制造教育部重點實驗室，山東濟(jì)南250061

0 引 言

在勢流理論下，海洋工程中的許多問題都可以歸結(jié)為邊界值問題，如波浪能發(fā)電裝置、海洋石油鉆井平臺，以及船舶與波浪、水流的相互作用等。在現(xiàn)有的求解方法中：采用解析法求解精確、計算速度快，但對求解域的幾何形狀要求較為苛刻，只適用于一些幾何形狀較為簡單的規(guī)則求解域，比如特征函數(shù)匹配展開法適用于截面為矩形的浮體與波浪的相互作用［1］，多極子法適用于截面為圓形的浮體與波浪的作用［2］；數(shù)值方法中的有限元法可以對復(fù)雜形狀進(jìn)行求解，但需要對整個求解域進(jìn)行離散，計算量較大［3-4］。

相比于需要對整個求解域進(jìn)行離散的方法，如有限元法、有限差分法等，邊界元法具有顯著的優(yōu)勢［5］。邊界元法最顯著的特點之一是在計算時更小的計算量和數(shù)據(jù)存儲；此外，邊界元法的數(shù)值精度通常要優(yōu)于有限元法。采用邊界元法僅需在求解域邊界上進(jìn)行離散，將三維問題轉(zhuǎn)化為二維問題，二維問題轉(zhuǎn)化為一維問題，即能很方便地處理無限域問題，這也是邊界元法在水動力學(xué)問題中能得到廣泛應(yīng)用的原因之一［6］。

本文將分別采用二維下的非連續(xù)參數(shù)化單元和參數(shù)化單元邊界元法解勢流速度場問題，以便在較低的單元數(shù)下得到較為理想的求解精度，并采用一種三維下的參數(shù)化單元邊界元法解勢流速度場問題，以便快速得到可以接受的平均誤差精度。

1 數(shù)學(xué)模型及積分方程建立

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

勢流問題中的流體為無旋、無粘性、不可壓縮的理想流體，流場域滿足拉普拉斯方程：

流場域滿足相應(yīng)的邊界條件如下：

式中：n為物面某點的法向向量，垂直于物面向外；?為速度勢；f為給定的函數(shù)，下標(biāo)1，2表示不同的給定函數(shù)（以下同）。

式（2）中的法向?qū)?shù)由式（3）定義：

式中，nx，ny分別為物面單位法向向量在x，y軸上的分量，法向量方向垂直物面向外。需注意的是，單位法向量在不同的分量上是不同的，其是關(guān)于x和y的函數(shù)。求解域的控制方程（1）和邊界條件已知，便構(gòu)成了邊界值問題，且存在特解。

1.2 積分方程的建立與離散

控制方程［7］存在基本解：

式中，ξ和η為選定點的坐標(biāo)。根據(jù)格林公式，容易得邊界積分方程

式中：S為曲線積分；C為積分路徑。

在邊界積分方程中，λ(ξ，η)的定義如下［8］：

將求解域的邊界L使用非連續(xù)參數(shù)化單元近似為L1+L2+…+LN，各離散單元以逆時針順序排列于求解域的邊界上。離散單元L1，L2，…，LN由求解域邊界上的點 (x1，y1)，(x2，y2)，…，(xN，yN)分隔，且N為有限值。各離散單元的邊界條件由式（2），得

每個單元的端點坐標(biāo)分別是(xk，yk)和(xk+1，yk+1)，且需注意，此處的 (xN+1，yN+1)=(x1，y1)。

貫穿整個非連續(xù)參數(shù)化單元，采用單元上點的f取值來近似fk。為了進(jìn)行線性近似，需要在單元上的2個不同點取fk的近似值。這里選取 2 個 點 (ξk，ηk) 和 (ξN+k，ηN+k) ，它 們 與 點(xk，yk)，(xk+1，yk+1)之間的長度均為τlk。其中，lk為單元k的長度，τ為一給定系數(shù)（0＜τ＜1/2）。單元k的具體結(jié)構(gòu)參見圖1。

圖1 非連續(xù)參數(shù)化單元示意圖Fig.1 Schematic diagram of discontinuous parameterized element

點 (ξk，ηk)，(ξN+k，ηN+k)處的?值分別由?k和?N+k指代。根據(jù) Ang［9］的研究，為了用?k和?N+k近似表示出單元上?值的線性變化，定義如下：

其中，(x，y)屬于Lk（Lk為第k個離散單元）。由式（8）可知，s(x，y)為Lk上的點 (x，y)與端點 (xk，yk)間 的 距 離 。 注 意 ，點 (ξk，ηk) ，(ξN+K，ηN+k) 的s(x，y)值分別為τlk和 (1-τ)lk。所需?(x，y)值的線性近似可以以?k，?N+k的形式給出，見式（9）。類似地，有，見式（10）。

當(dāng) 0＜τ＜1/2時，近似表達(dá)式（9）、式（10）定義為非連續(xù)單元。在對求解域邊界L進(jìn)行離散時，任意元素都有?k和pk，且其中有且僅有一個已知，如邊界上?k已知，pk則未知。因此，式（9）、式（10）中有 2N個未知項。將式（9）和式（10）代入式（5），可得

其中，

式中，Ck為單元積分路徑。

為了求解式（11）右側(cè)的2N個未知項，建立了一個2N階線性代數(shù)方程組，即可以依次取點(ξ，η)為 (ξm，ηm)，其中m=1，2，…，2N，則又有式（16）。由圖 1，可知點 (ξk，ηk)和 (ξN+k，ηN+k)為單元內(nèi)的兩點，則λ(ξm，ηm)=1/2，m=1，2，…，2N。式（16）又可寫為式（17），式中的zk，zN+k表示元素上的未知量。其中在邊界單元上，當(dāng)?k已知時，amk，am(N+k)，bmk分別由式（18）、式（19）和式（20）給出，δmk為Kronecker數(shù)。在邊界單元上，當(dāng)pk已知時，情況與上述類似，此處不再贅述。(其中p=1，2，3，4)系數(shù)的參數(shù)化計算方法參見文獻(xiàn)［9］。

其中，

2 二維勢流問題算例

2.1 矩形勢流場算例

第1個算例，考慮矩形勢流場問題，定解問題如下：

由文獻(xiàn)［10］可知，上述問題的解析解為

如圖2所示，在分析計算中，將矩形計算域邊界離散為30個單元，將邊界積分方程離散，可以得到一個線性方程組，解此方程組即可求得所需的速度勢或其他未知量。

圖2 矩形流場示意圖Fig.2 Schematic diagram of square flow field

在劃分30個單元、采用非連續(xù)參數(shù)化單元離散的情況下，速度勢的最大相對誤差為0.63%，平均相對誤差為0.19%，具體計算結(jié)果如表1所示；在劃分單元數(shù)加倍、采用連續(xù)常數(shù)單元的情況下，速度勢的最大相對誤差為70.86%，平均相對誤差為8.42%，計算結(jié)果如表2所示。

2.2 圓形勢流場算例

第2個算例，考慮圓形勢流場問題，定解問題如下：

式中:r，θ為極坐標(biāo)中的變量；r0為圓形流場域的半徑，此處r0=1。其解析解為：

如圖3所示，在分析計算中，將圓形計算域邊界離散為32個單元，將邊界積分方程離散，亦可得到一個線性方程組，解此方程組即可求得所需的速度勢或其他未知量。

表2 參數(shù)化單元速度勢計算結(jié)果（案例1）Table 2 The calculation results of velocity potential by parameterized elements（case 1）

在劃分32個單元、采用非連續(xù)參數(shù)化單元離散的情況下，速度勢的最大相對誤差為0.02%，平均相對誤差為0.01%，計算結(jié)果如表3所示；在劃分單元數(shù)加倍、采用連續(xù)常數(shù)單元的情況下，速度勢的最大相對誤差為0.97%，平均相對誤差為0.08%，計算結(jié)果如表4所示。

圖3 圓形流場示意圖Fig.3 Schematic diagram of circle flow field

2.3 非連續(xù)參數(shù)化單元與參數(shù)化單元比較

在二維情況下，分別采用非連續(xù)參數(shù)化單元和參數(shù)化單元邊界元法解勢流速度場問題，其中非連續(xù)參數(shù)化單元可以在較低的單元數(shù)下得到較為理想的求解精度。

由于非連續(xù)參數(shù)化單元上的物理量并非與單元節(jié)點處近似，故相比于參數(shù)化單元邊界元法，非連續(xù)參數(shù)化單元在求解多邊界問題時不但具有常數(shù)單元交界點無需特殊處理的優(yōu)點，而且還具有線性及高次單元在較少單元數(shù)下得到較高求解精度的特點。

當(dāng)采用傳統(tǒng)常數(shù)或線性單元的直接邊界元法解勢流速度場問題時，除了邊界近似帶來的誤差，還有在求矩陣系數(shù)時由高斯積分法近似引起的誤差。非連續(xù)參數(shù)化單元由于是將單元上的點進(jìn)行參數(shù)化，求解系數(shù)時使用解析表達(dá)式求解，即柯西主值積分，故其只有邊界近似誤差以及系統(tǒng)矩陣求解時帶來的誤差。

表3 非連續(xù)參數(shù)化單元速度勢計算結(jié)果（案例2）Table 3 The calculation results of velocity potential by discontinuous parameterized elements（case 2）

表4 參數(shù)化單元速度勢計算結(jié)果（案例2）Table 4 The calculation results of velocity potential by parameterized elements（case 2）

3 三維勢流問題算例

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

考慮無窮遠(yuǎn)邊界的流場［11-12］，有沿x軸負(fù)向的、速度為1的均勻來流。流體無旋、無粘性、不可壓縮，流場中有一個半徑為1的圓球。流場速度勢由下式表示［6］：

式中：Φ0為總速度勢；V∞為無窮遠(yuǎn)處的來流速度；φ為圓球的擾動速度勢。將求解的直接變量選擇為擾動速度勢φ，擾動速度勢φ滿足定解方程

式中：Ω為流場區(qū)域；S為物面邊界；nQ為Q點的單位法向量。

將基本解?(x，y)=-代入三維邊界積分方程，離散并建立方程組，求解方程組后，便可得單元上的未知量。

3.2 離散與數(shù)值求解

使用三角形常數(shù)單元對求解域表面進(jìn)行離散，并使用一種對單元上的點進(jìn)行參數(shù)化的方法對問題進(jìn)行求解。邊界積分方程離散后，其求解過程與二維非連續(xù)參數(shù)化單元類似，此處不再贅述。圖4為圓球劃分三角網(wǎng)格效果圖。

在這里，使用式（28）對單元上的點進(jìn)行參數(shù)化：

式中：Xk，Yk，Zk為第k個單元上角點坐標(biāo)的參數(shù)化表示；u，v為變換坐標(biāo)中的量。

式中，(x1，y1，z1)，(x2，y2，z2)，(x3，y3，z3)分別為三角形元素上角點1，2，3的坐標(biāo)。

式中：Φ3D為三維拉普拉斯方程的基本解；Jk為雅可比常數(shù)，由下式給出：

其中，

式（32）和式（33）又可以寫為式（35）及式（36），其中的t為前面所述的u。

這樣，就可以使用高斯積分公式（37）對上式進(jìn)行數(shù)值求解［9］。

圖4 圓球劃分三角網(wǎng)格效果Fig.4 The effect of sphere divided into triangular grids

3.3 速度求解

求得物面控制點勢函數(shù)后，使用文獻(xiàn)［6］所采用的方法計算物面的速度分布，詳細(xì)步驟參見文獻(xiàn)［13］。

本文選取不同的單元數(shù)，分別采用上述方法進(jìn)行了數(shù)值實驗，以參數(shù)化單元高斯積分法計算系數(shù)矩陣的結(jié)果如圖5～圖10所示。圖中，相對誤差即為相對解析解［11-12］的誤差。

當(dāng)單元數(shù)為4 512時，速度勢的最大誤差為19.76%，平均誤差為1.72%；速度的最大誤差為29.78%，平均誤差為5.18%。

圖5 4 512個網(wǎng)格速度勢計算結(jié)果Fig.5 The calculation results of velocity potential of 4 512 grids

圖6 4 512個網(wǎng)格速度計算結(jié)果Fig.6 The calculation results of velocity of 4 512 grids

圖7 4 512個網(wǎng)格速度矢量結(jié)果Fig.7 The velocity vector results of 4 512 grids

圖8 9 112個網(wǎng)格速度勢計算結(jié)果Fig.8 The calculation results of velocity potential of 9 112 grids

圖9 9 112個網(wǎng)格速度計算結(jié)果Fig.9 The calculation results of velocity of 9 112 grids

基于網(wǎng)格的劃分方法，球的兩個極點附近的單元數(shù)較少，即便增加總網(wǎng)格數(shù)，誤差也較其他位置節(jié)點處的突出。由于所使用高斯積分公式的數(shù)值積分精度不高，因此，三維問題下的計算結(jié)果不如二維問題下的好，可以考慮結(jié)合采用改進(jìn)網(wǎng)格劃分和采用非連續(xù)單元計算，也可以采用其他數(shù)值積分方法來提高計算精度。

圖10 9 112個網(wǎng)格速度矢量結(jié)果Fig.10 The velocity vector results of 9 112 grids

速度幅值計算結(jié)果的精度是可以接受的，但該計算結(jié)果也未能精確反映解析解給出的速度變化趨勢。由圖可以看到，速度勢和速度的平均相對誤差都可以接受，但某些節(jié)點的誤差偏大。這主要是由引入的高斯積分公式精度不高、數(shù)值計算誤差所致；另外，不僅速度勢函數(shù)計算本身存在誤差，物面速度計算也存在誤差，故可以看到速度勢與速度的誤差特性有一定的相關(guān)性。

當(dāng)單元數(shù)為9 112個時，速度勢最大誤差為17.62%，平均誤差為1.23%；速度最大誤差為27.98%，平均誤差為4.92%。

數(shù)值實驗表明，相對于文獻(xiàn)［6］的7點高斯積分法，本文方法的優(yōu)點是計算同數(shù)量網(wǎng)格所需時間更少，速度更快；缺點是求解精度存在問題，平均精度尚在可接受范圍內(nèi)，但某些節(jié)點的計算誤差偏大。通過文中給出的2種情況可以看出，在網(wǎng)格加密的情況下，本文方法的精度提高并不大，數(shù)值實驗表明，在約4 000網(wǎng)格數(shù)時就已達(dá)到相當(dāng)于10 000網(wǎng)格數(shù)時的精度，也就是說，求解精度在 4 000網(wǎng)格數(shù)左右達(dá)到最優(yōu)結(jié)果；在計算機(jī)可計算范圍內(nèi)，隨著網(wǎng)格的不斷增加，計算結(jié)果中會出現(xiàn)誤差非常大的異常點，筆者認(rèn)為，這主要是由引入高斯積分公式及參數(shù)化單元計算時進(jìn)行的坐標(biāo)變換引起單元節(jié)點與實際節(jié)點不一致所引起。

4 結(jié) 語

邊界元法的計算量主要集中在影響系數(shù)矩陣的計算上，誤差也大多源于此。本文對多種參數(shù)化單元方法進(jìn)行了數(shù)值實驗，分析了各種方法對最終結(jié)果誤差的影響。

二維下的非連續(xù)參數(shù)化單元和參數(shù)化單元邊界元法，由于單元上的節(jié)點參數(shù)化，矩陣系數(shù)計算均為解析式計算，即柯西主值積分，速度快、精度高，沒有引起數(shù)值計算誤差。其中，非連續(xù)參數(shù)化單元在求解多邊界問題時不但兼顧了常數(shù)單元交界點無需特殊處理，而且兼顧了二階及以上單元在較少單元數(shù)下得到較高求解精度的優(yōu)點。

三維下的參數(shù)化單元邊界元法，由于單元上的節(jié)點參數(shù)化，計算速度相比于其他方法稍快，但矩陣系數(shù)計算引入了高斯積分公式，引進(jìn)了數(shù)值計算誤差，在計算系數(shù)積分時進(jìn)行了坐標(biāo)變換，使得變換后求得的單元節(jié)點與實際節(jié)點存在偏差，故又引入了新的誤差。數(shù)值實驗表明，計算平均誤差可以接受，但存在著某些節(jié)點的誤差偏大的問題，易導(dǎo)致求解不夠精確?？梢酝ㄟ^改進(jìn)網(wǎng)格劃分和采用非連續(xù)單元計算，或其他可通過坐標(biāo)轉(zhuǎn)換進(jìn)行解析計算的單元來計算。

本文的二維非連續(xù)參數(shù)化單元可以用于對波浪與浮體，或者潛體的作用機(jī)理進(jìn)行研究，單元數(shù)少，求解精度高；三維參數(shù)化單元雖然求解精度差，但求解速度快，若進(jìn)一步改進(jìn)，如非連續(xù)單元在三維時，可提高運(yùn)算精度。

本文所涉及的參數(shù)化單元邊界元法的編程與傳統(tǒng)邊界元法的編程程序相比具有通用性，可以解決一類在勢流框架下的工程及科學(xué)問題。