王小軍，劉曉宇，杜亞楠

(重慶大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，重慶 400030)

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

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

1 復(fù)連通區(qū)域Stokes問題的邊界積分方程

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

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

圖1 研究區(qū)域

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

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

2 變分方程的形成及處理

2.1 形成變分方程

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

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

2.2 離散變分方程

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

其中:

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

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

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

3 程序設(shè)計

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

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

在輸入功能子程序INPUT中要實現(xiàn)的功能是輸入邊界節(jié)點坐標(biāo)、已知邊界條件和所需計算的內(nèi)點坐標(biāo)。

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

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

計算4點Gauss數(shù)值積分公式的子程序，調(diào)用為Gauss數(shù)值積分作準(zhǔn)備的子程序，將坐標(biāo)轉(zhuǎn)變?yōu)闊o因次積分點坐標(biāo)，實現(xiàn)非奇異情況下的數(shù)值積分。

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

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

輸出結(jié)果的子程序?qū)⒂嬎愠龅慕Y(jié)果寫入文本文件，以便對比和畫流線圖。

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

4 數(shù)值算例

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

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

圖2 算例1流線圖

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

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

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

圖3 算例2流線圖

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

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

圖4 算例3的漩渦圖

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

［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］.重慶大學(xué)學(xué)報，2006(2):29－30.

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

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

［5］祝家麟.定常Stokes問題的邊界積分方程法［J］.計算數(shù)學(xué)，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.