王小軍，劉曉宇，杜亞楠

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

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

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

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

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

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

圖1 研究區(qū)域

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

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

2 變分方程的形成及處理

2.1 形成變分方程

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

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

2.2 離散變分方程

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

其中:

把上述線(xiàn)性方程組寫(xiě)成簡(jiǎn)潔形式:AU=B，其中:

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

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

3 程序設(shè)計(jì)

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

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

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

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

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

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

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

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

輸出結(jié)果的子程序?qū)⒂?jì)算出的結(jié)果寫(xiě)入文本文件，以便對(duì)比和畫(huà)流線(xiàn)圖。

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

4 數(shù)值算例

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

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

圖2 算例1流線(xiàn)圖

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

計(jì)算矩形內(nèi)單位圓以外的點(diǎn)的流速并畫(huà)出流線(xiàn)圖。

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

圖3 算例2流線(xiàn)圖

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

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

圖4 算例3的漩渦圖

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

［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é)報(bào)，2006(2):29－30.

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

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

［5］祝家麟.定常Stokes問(wèn)題的邊界積分方程法［J］.計(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.