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

陳 瑤，郭元輝

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

0 引言

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

Δφ=0.

(1)

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

1 問(wèn)題及方程

圖1 示意圖Fig.1 The diagram

假設(shè)在圖1所示的兩個(gè)平行平板間，通入速度為2 m/s的空氣，出口敞開(kāi).平行平板的長(zhǎng)度為10 m，寬度為5 m.試求流動(dòng)區(qū)域內(nèi)的速度式分布[7].

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

(2)

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

(3)

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

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

2 有限元理論

定義1.2對(duì)于非負(fù)整數(shù)k和實(shí)數(shù)p≥1,定義

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

這是一個(gè)Banach空間.賦予范數(shù)

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

(4)

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

(5)

(6)

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

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

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

連續(xù)性|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.

定理的證明用能量范數(shù)和Ritz表現(xiàn)定理可以得到，此處不加證明地引用.

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

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

ii.P是K上的有限維函數(shù)空間；

iii.N={N1,...,Nn}是一組由節(jié)點(diǎn)P'構(gòu)成的基.

(9)

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

證明：i. 對(duì)任意g∈L2(Ω),問(wèn)題

(10)

ii.存在常數(shù)C滿足

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

(11)

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

=(u-uh,g)

=a(u-uh,φg)

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

證畢.

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

根據(jù)有限元理論，通過(guò)程序設(shè)計(jì)實(shí)現(xiàn)Laplace方程的計(jì)算，這里有幾個(gè)主要步驟.

i.區(qū)域的剖分.通常區(qū)域可以使用三角形剖分或四邊形剖分，這里根據(jù)題目條件，選擇四邊形剖分，給出若干行列，把區(qū)域劃分成大小相等的小四邊形.標(biāo)記上節(jié)點(diǎn)和單元編號(hào).

ii.選擇基函數(shù).因?yàn)閷?duì)每一個(gè)單元，只需要用到四個(gè)節(jié)點(diǎn)的函數(shù)值，所以插值函數(shù)用雙線性形式：

φ(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.計(jì)算結(jié)果.u=A-1B.

通過(guò)Matlab編程，計(jì)算結(jié)果如表1所示.

表1 區(qū)域剖分為5×10個(gè)四邊形時(shí)，節(jié)點(diǎn)所對(duì)應(yīng)的uTab.1 u corresponding to the node when the region is divided into50quadrilaterals

4 使用PDE-Tool與freeFem++計(jì)算

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

同樣，為了驗(yàn)證計(jì)算的結(jié)果，這里也可以使用網(wǎng)上的免費(fèi)軟件Freefem++[14]，它是一個(gè)高度集成化的有限元軟件，其編程簡(jiǎn)單、直觀、高效，其核心部分的書(shū)寫(xiě)與變分形式一一對(duì)應(yīng)，并且把第一邊界、第二邊界的處理也封裝了.例如：

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++的計(jì)算結(jié)果如圖3所示.

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