復連通二維Stokes問題的Galerkin邊界元解法

王小軍，劉曉宇，杜亞楠

(重慶大學數學與統計學院，重慶 400030)

邊界元方法是求解不可壓縮粘性流體Stokes流動問題數值解的理想方法。對于二維問題，在區域的邊界是封閉曲線的情況下，如無界流體的繞流問題或有界固壁容器中流體的流動問題等Stokes問題，用邊界元方法求解已有豐富的研究成果，然而多側重于理論分析。本文擬對復連通二維Stokes問題的Galerkin邊界元解法中的程序設計等技術性問題進行探討。

根據Green公式和Stokes算子的基本解可以推導出對應于復連通閉曲線Stokes問題的邊界積分方程，得出基于單層位勢的間接邊界積分方程，這是一個第1類的Fredholm向量積分方程。對與之等價的邊界變分方程采用Galerkin邊界元求解，得出單層位勢的向量密度，進而得出流場中任意點的流速值。本文主要討論其中用利用Fortran計算所求點的流速與用Matlab畫出相應的流線圖的主要過程。

1 復連通區域Stokes問題的邊界積分方程

設Γ1是R2中有限區域Ω的封閉曲線，Γ2是Γ1形成的封閉區域中的一條封閉曲線，Γ2形成的區域為 Ω0+，Ω0－=Ω －Ω0+。設 Γ=Γ1+Γ2，考慮如下的問題:

其中:u=( u1，u2)是流速;p是壓強;u0是閉曲線邊界Γ1以及內部閉曲線邊界Γ2上的已知函數。

圖1 研究區域

設想包含在Ω內的不可壓縮粘性流體整體嵌入在平面無限域上的不可壓縮粘性流體之中。本文只計算有界域上的流動，故假定在無窮遠處流速為零。由參考文獻［5］，經計算可得到基于單層位勢的邊界積分方程:

其中 t= ( t1，t2)是待定的密度函數，其分量ti是穿越Γ時的躍變。類似地可推出流體壓力場的積分表達式，但本文約去壓力場的計算。

2 變分方程的形成及處理

2.1 形成變分方程

由單層位勢出發來研究邊界量u0和邊界量t的關系，由于u在穿越邊界時的連續性，從式(2)得到聯系已知函數 u0和中間變量 t的積分方程組:

對于方程(3)，方程兩邊同乘以 t '(y)，然后在Γ上積分，便得到變分方程組:

2.2 離散變分方程

為了數值求解變分方程(4)，邊界Γ必須離散為一系列單元，從而把邊界積分方程轉化為線性代數方程組進行求解。設把邊界Γ1離散為n1個單元，有單元端點n1個;把邊界Γ2離散為n2個單元，有單元端點 n2個，則邊界 Γ共劃分為 n( n =n1+n2)個單元，共有n個單元端點。在實踐中，邊界單元一般被離散為常單元、線性單元或者高次單元，也有采用樣條插值方式的邊界單元。在本文中，采用線性單元就二維問題進行計算。基函數選取如下:

其中:

把上述線性方程組寫成簡潔形式:AU=B，其中:

2.3 變分方程中重積分的處理

要求解上述線性方程組(5)，最主要的問題是求解系數矩陣。然而由于積分存在奇異性，若采用數值積分，必然產生較大的誤差。為了避免出現這樣的問題，在求二重積分時，第1重積分存在奇異性時采用精確解析積分，不存在奇異性時以及第2重積分采用Gauss數值積分。在利用Gauss數值積分計算第2重積分時，對于以閉曲線為邊界的區域問題而言，基函數無奇異性，可以方便的使用通常的Gauss數值積分公式。解析公式的推導方法與具體公式見參考文獻［2］。利用邊界積分方程得到密度函數后，將其數值帶入解的表達式即得流速。

3 程序設計

鑒于不同程序設計語言的優勢，程序采用兩種語言，首先用Fortran計算出所求點的流速，保存結果在文本文件中，然后用Matlab讀取畫出流線圖。

Fortran程序采用模塊化設計，包括主程序MIAN以及8個子程序，其中主程序控制程序的執行順序，采用多個輸入輸出文本文件讀取寫入以便于與Matlab結合。在主程序中設置全局變量，包括了邊界點，邊界點對應的流速，邊界單元長度，稀疏矩陣，右端項，單層位勢的向量密度，所求點的坐標以及所求解。

在輸入功能子程序INPUT中要實現的功能是輸入邊界節點坐標、已知邊界條件和所需計算的內點坐標。

在形成系數矩陣的子程序FMAT中，首先通過兩層循環采用4點Gauss數值積分公式形成右端項，然后通過3層循環調用解析積分公式和4點Gauss數值積分公式形成系數矩陣，其中在第1層積分中，存在奇異性時調用解析積分公式，不存在奇異性時調用4點Gauss數值積分公式。

計算解析積分公式的子程序FUN中，通過判斷語句根據積分源點到有向線段的距離是否為零分情況計算系數矩陣對應元素的數值。

計算4點Gauss數值積分公式的子程序，調用為Gauss數值積分作準備的子程序，將坐標轉變為無因次積分點坐標，實現非奇異情況下的數值積分。

計算線性代數方程組的源程序，采用列主元消去法求解邊界積分方程組。

計算流速的子程序中采用2層循環，將求得的密度數值帶入解的表達式求得未知點的流速及其方向。

輸出結果的子程序將計算出的結果寫入文本文件，以便對比和畫流線圖。

Matlab程序中首先對區域進行剖分，得到所求點坐標，將坐標輸出到文本文件，由Fortran程序計算出流速后，加載并讀取保存結果的文本文件，然后調用已知函數畫出流線圖，最后畫出邊界并確定坐標名稱。

4 數值算例

算例1 單位圓流速為(0，0)，矩形區域［－5，5，－10，10］邊界上的速度為(1，0)，計算矩形內單位圓以外的點的流速并畫出流線圖。

在此算例中，內部圓邊界劃分為16個單元，外部閉曲線劃分為60個單元，邊界總剖分為76單元，由Matlab剖分在矩形區域內確定 1 521個點，除去圓內部余1 492個有效點，由Fortran程序計算出流速后輸出到文本文件，耗費時間以秒為單位，再由Matlab程序讀取畫出流線圖，耗費時間也是以秒為單位，與其他方法計算結果相對比精度較高，畫出的流線圖與實際情況符合很好(圖2)。

圖2 算例1流線圖

算例2 單位圓流速為(0，0)，矩形區域［－5，5，－10，10］邊界上的速度為

計算矩形內單位圓以外的點的流速并畫出流線圖。

在此算例中，邊界剖分與計算流程與算例1完全一致，計算、畫圖耗用時間以及結果精確度也完全一致(圖3)。

圖3 算例2流線圖

算例3 單位圓流速為(0，0)，矩形區域［－5，5，－10，10］邊界上的速度為計算矩形內單位圓以外的點的流速并畫出流線圖。

在此算例中，邊界剖分與計算流程與上述兩例完全一致，計算、畫圖耗用時間以及結果精確度也完全一致。不同點在于，Matlab畫圖中放大后將清晰的出現2個小的漩渦。漩渦如圖4所示。

圖4 算例3的漩渦圖

通過以上數值算例可以看出，該種方法精度較高，計算結果符合客觀事實。

［1］Hsiao G C.A modified Galerkin scheme for equations with natural boundary conditions［C］//Vichnevetsky R，Vigness J.Numerical Mathematics and Applications.B.V:Elsevier Science Publishers，1986:193 －197.

［2］向瑞銀.平面定常Stokes方程的Galerkin邊界元解法［J］.重慶大學學報，2006(2):29－30.

［3］張耀明，溫衛東.平面定常Stokes問題的無奇異第一類邊界積分方程［J］.計算數學，2005，2(1):1－10.

［4］祝家麟.橢圓邊值問題的邊界元分析［M］.北京:科學出版社，1991.

［5］祝家麟.定常Stokes問題的邊界積分方程法［J］.計算數學，1986(8):281－289.

［6］Zhu Jialin.The boundary integral equation method for solving stationary Stokes problems［C］//Feng Kang.Proc.of the 1984 Beijing Symp.On Diff.Geometry and Diff.Equations.Beijing:Science Press，1985.