參數化單元邊界元法解勢流速度場問題

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

1山東大學海洋研究院，山東濟南250100

2山東大學機械工程學院，山東濟南250061

3山東大學高效清潔機械制造教育部重點實驗室，山東濟南250061

0 引 言

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

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

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

1 數學模型及積分方程建立

1.1 數學模型

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

流場域滿足相應的邊界條件如下：

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

式（2）中的法向導數由式（3）定義：

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

1.2 積分方程的建立與離散

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

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

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

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

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

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

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

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

點 (ξk，ηk)，(ξN+k，ηN+k)處的?值分別由?k和?N+k指代。根據 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）。

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

其中，

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

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

其中，

2 二維勢流問題算例

2.1 矩形勢流場算例

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

由文獻［10］可知，上述問題的解析解為

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

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

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

2.2 圓形勢流場算例

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

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

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

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

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

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

2.3 非連續參數化單元與參數化單元比較

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

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

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

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

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

3 三維勢流問題算例

3.1 數學模型

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

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

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

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

3.2 離散與數值求解

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

在這里，使用式（28）對單元上的點進行參數化：

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

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

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

其中，

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

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

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

3.3 速度求解

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

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

當單元數為4 512時，速度勢的最大誤差為19.76%，平均誤差為1.72%；速度的最大誤差為29.78%，平均誤差為5.18%。

圖5 4 512個網格速度勢計算結果Fig.5 The calculation results of velocity potential of 4 512 grids

圖6 4 512個網格速度計算結果Fig.6 The calculation results of velocity of 4 512 grids

圖7 4 512個網格速度矢量結果Fig.7 The velocity vector results of 4 512 grids

圖8 9 112個網格速度勢計算結果Fig.8 The calculation results of velocity potential of 9 112 grids

圖9 9 112個網格速度計算結果Fig.9 The calculation results of velocity of 9 112 grids

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

圖10 9 112個網格速度矢量結果Fig.10 The velocity vector results of 9 112 grids

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

當單元數為9 112個時，速度勢最大誤差為17.62%，平均誤差為1.23%；速度最大誤差為27.98%，平均誤差為4.92%。

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

4 結 語

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

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

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

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

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