陳 瑤，郭元輝

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

0 引言

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

Δφ=0.

(1)

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

1 問題及方程

圖1 示意圖Fig.1 The diagram

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

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

(2)

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

(3)

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

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

2 有限元理論

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

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

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

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

(4)

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

(5)

(6)

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

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

定理1.1[5](Lax-Milgram Lemma)假設(shè)V是一個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有限元是一個具有如下性質(zhì)的三元組(K,P,N)

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

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

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

(9)

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

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

(10)

ii.存在常數(shù)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 程序設(shè)計方法

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

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

ii.選擇基函數(shù).因為對每一個單元，只需要用到四個節(jié)點的函數(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.計算結(jié)果.u=A-1B.

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

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

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

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

同樣，為了驗證計算的結(jié)果，這里也可以使用網(wǎng)上的免費軟件Freefem++[14]，它是一個高度集成化的有限元軟件，其編程簡單、直觀、高效，其核心部分的書寫與變分形式一一對應(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++的計算結(jié)果如圖3所示.

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