理想流體的Laplace方程有限元解法

陳 瑤，郭元輝

(1.西華師范大學數學與信息學院，四川南充 637009；2.西華師范大學教育信息技術中心，四川南充 637009)

0 引言

在大自然現象中，科學的研究，往往借助于各種復雜的偏微分方程的描述.拉普拉斯方程是最常見的一種，它廣泛應用于流體力學、彈性力學、熱傳導，電磁波、現代光學等.它的基本形式為[1,2]

Δφ=0.

(1)

如何求解偏微分方程是現代科學計算所研究的一個主要問題.科學計算與實驗、理論已并列成為現代科學研究的三種方法.科學計算不僅是一種數字計算的手段，而且也是一種研究方法[3].從早期的古典解、解析解，到近代的數值解、近似解，有限差分、有限元[4-6]、邊界元及深度學習、人工智能算法等，各種方法不斷出現、創新，帶來計算科學和技術的發展，以及更高的數值精度和數學理論、科學理論的進步.

1 問題及方程

圖1 示意圖Fig.1 The diagram

假設在圖1所示的兩個平行平板間，通入速度為2 m/s的空氣，出口敞開.平行平板的長度為10 m，寬度為5 m.試求流動區域內的速度式分布[7].

理想流體[8]通常定義為其內部沒有摩擦的流體，稱為無黏性流體.這里，假設空氣滿足理想流體的基本特性，也就是黏滯系數η為0；同時，不可壓縮，密度ρ為常數[9,10].那么，理想流體的連續性方程為

(2)

其中，u,v分別表示速度勢函數φ(x,y)在x,y方向的導數?φ/?x,?φ/?y，代入式(2)得到

(3)

這就是Laplace方程形式.本題中有兩類邊界條件.

第一類邊界，也叫Dirichlet邊界條件，在邊界上知道速度勢φ(x,y)的值.第二類邊界，又叫Neumann邊界條件，知道邊界上速度勢的法向導數?φ/?n.

2 有限元理論

定義1.2對于非負整數k和實數p≥1,定義

Wk,p(Ω)={u∈Lp(Ω):??u∈Lp(Ω),?|α|k}.

這是一個Banach空間.賦予范數

設f?L2(Ω)，這里考慮更一般的Poisson方程，本題Laplace方程僅僅是f=0的特例.

(4)

設u∈H2是方程(4)的古典解，定義空間V:={v∈H1(Ω):v|Γ=0}，

(5)

(6)

這樣求出的弱解與古典解具有等價性[4].v所在的空間，稱為試驗空間；u所在的空間成為容許空間.當使用Galerkin方法時，兩者取同一空間，稱為能量空間.

那么，如何保證弱解的存在唯一性？

定理1.1[5](Lax-Milgram Lemma)假設V是一個Hilbert空間，定義范數‖ . ‖和內積(. , .)，a(. , .)是映射V×V→R上的雙線性函數，存在常數α，β>0滿足

連續性|a(u,v)|β‖u‖‖v‖,?u,v∈V,

(7)

和V-橢圓a(v,v)≥α‖v‖2,?v∈V.

(8)

那么，存在唯一的解u∈V，滿足a(u,v)=F(v),?v∈V.

定理的證明用能量范數和Ritz表現定理可以得到，此處不加證明地引用.

定義1.3有限元是一個具有如下性質的三元組(K,P,N)

i.K∈d是一個具有分片光滑邊界的閉集；

ii.P是K上的有限維函數空間；

iii.N={N1,...,Nn}是一組由節點P'構成的基.

(9)

‖u-uh‖L2(Ω)Ch2|u|H2(Ω).

證明：i. 對任意g∈L2(Ω),問題

(10)

ii.存在常數C滿足

‖φg‖H2(Ω)C‖g‖L2(Ω).

(11)

iii.取g=u-uh,讓φg是方程(9)的解，那么

=(u-uh,g)

=a(u-uh,φg)

=a(u-uh,φg-Ihφg)

證畢.

3 程序設計方法

根據有限元理論，通過程序設計實現Laplace方程的計算，這里有幾個主要步驟.

i.區域的剖分.通常區域可以使用三角形剖分或四邊形剖分，這里根據題目條件，選擇四邊形剖分，給出若干行列，把區域劃分成大小相等的小四邊形.標記上節點和單元編號.

ii.選擇基函數.因為對每一個單元，只需要用到四個節點的函數值，所以插值函數用雙線性形式：

φ(x,y)=ax+by+cxy+d.

for k=1:dys

dybh=dyjd(k,:); dyzb=jdzb(dybh,:);

xc=sum(dyzb(:,1))/4;yc=sum(dyzb(:,2))/4;

b=(dyzb(2,1)-dyzb(1,1))/2;c=(dyzb(4,2)-dyzb(1,2))/2;

ke(1,1)=(b^2+c^2)/(3*b*c);ke(1,2)=(b^2-2*c^2)/(6*b*c);

ke(1,3)=-(b^2+c^2)/(6*b*c);ke(1,4)=(c^2-2*b^2)/(6*b*c);

ke(2,1)=ke(1,2);ke(2,2)=ke(1,1);ke(2,3)=ke(1,4);ke(2,4)=ke(1,3);

ke(3,1)=ke(1,3);ke(3,2)=ke(2,3);ke(3,3)=ke(1,1);ke(3,4)=ke(1,2);

ke(4,1)=ke(1,4);ke(4,2)=ke(2,4);ke(4,3)=ke(3,4);ke(4,4)=ke(1,1);

kk(dybh,dybh)=kk(dybh,dybh)+ke; % 總裝

end

v.計算結果.u=A-1B.

通過Matlab編程，計算結果如表1所示.

表1 區域剖分為5×10個四邊形時，節點所對應的uTab.1 u corresponding to the node when the region is divided into50quadrilaterals

4 使用PDE-Tool與freeFem++計算

在MATLAB中，利用偏微分方程工具PDE-Tool，直接代入方程參數與邊界條件[13]，可以得到與表-1一樣的結果，如圖2所示.

同樣，為了驗證計算的結果，這里也可以使用網上的免費軟件Freefem++[14]，它是一個高度集成化的有限元軟件，其編程簡單、直觀、高效，其核心部分的書寫與變分形式一一對應，并且把第一邊界、第二邊界的處理也封裝了.例如：

problem laplace(u,v)=int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))

-int2d(Th)(Fh*v)

+on(2,u=0)

-int1d(Th,1)(0*v)-int1d(Th,3)(0*v)

-int1d(Th,4)(2*v);

使用Freefem++的計算結果如圖3所示.

圖2 PDE-Tool效果圖Fig.2 Calculation effect of PDE tool圖3 Freefem++計算效果圖Fig.3 Freefem++ calculation effect picture