二維三角形線性插值有限元法的原理與應(yīng)用

1、二維三角形單元的有限元程序編寫及應(yīng)用葛青宇1182201150 研電1803班一、二維有限元理論介紹1、電磁場邊值問題泊松方程是電磁場位函數(shù)電磁場邊值問題一般滿足不同的微分方程和相應(yīng)的邊界條件。 滿足方程中的一大類，如下所示：(2)W 2u = f.1 =U。；u a :n=bu c-2,3其中：(1)式稱為泊松方程， 為電磁場的位函數(shù)， 為表征空間介質(zhì)材料的常數(shù)，表示產(chǎn)生電磁場位函數(shù)的源；(2)式表示第一類邊界條件，為邊界上已知的位函數(shù)；(3)式為第三類邊界條件， 和 為已知的邊界條件系數(shù)。2、二維伽遼金加權(quán)余量方程以插值方式構(gòu)造近似函數(shù)，的近似解表示為nnU f M nun(4)n 1 ,

2、 Mnn ,相應(yīng)待定常數(shù)為U1,U2, ,Unn根據(jù)基函數(shù)和權(quán)函數(shù)與單元形狀函數(shù)的關(guān)系,在有限元法中一個單元內(nèi)基函數(shù)和權(quán)函數(shù)的值等于相應(yīng)的形狀函數(shù)的值，表不為在單兀e中，M n =N, M m =Nj。e Je, j采用伽遼金加權(quán)余量形式， 設(shè)一組與基函數(shù)相同的權(quán)函數(shù)Mi,M2,，Mnn ,將權(quán)函 數(shù)代入伽遼金加權(quán)余量方程，得一組方程 TOC o 1-5 h z -.Mm(au)d,J = Mmfd其中，m=1, 2,，nn,對應(yīng)節(jié)點的總體編號。對(4)式應(yīng)用格林定理，并將第二、 三類邊界條件代入得aMm； ud Mm(bu c)d- = Mmfd 1(6)QTQ其中m同上。將不含未知函數(shù)u的

3、項移至方程右端，得電磁場數(shù)值計算結(jié)課大作業(yè)a、Mm & udJ - QbMmud： = Mmfdj cMmd： 將由基函數(shù)線性組合構(gòu)成的近似函數(shù)代入方程，得.nnnn TOC o 1-5 h z atMm(MnUn)dJ-RbMmV (MnUjd： = (Mm f 2J (CMm)d：(8)1n3n鼻MJ將近似函數(shù)的梯度運算表示出來，得nnnnaVMmT C Mn)UndbMm (MnUn)d(9)(10)；-1n 1n 1=(Mmf)dn (cMm)d權(quán)函數(shù)與求和無關(guān)，可以將其移到求和號之內(nèi)，得nnnnaGMmJMn)Und-b (MmMn)Und： n 1n 4=(Mmf)d? (cMm)

4、d：將未知節(jié)點變量從積分中提出得：nnnn_-(aMm ； Mn)d-UneMmMn)d：Un = .(Mmf)dJ (cMm)d1(11)n“2 L可以看出代數(shù)方程如下所示：Km,lUi - Km,2U2+&,冊/0=%(12)對比式(10),可得：(13)(14)(15)(16)Km,n=.(alMmMn)dC-D(bMmMn)dRm = . (Mmf)d，(cMm)d：將(12)、(13)式化作單元上積分之和得形式，得：nenebKm,n =(aMm JMn2.(bMmMn)d1 eEeT EbnenebRm(Mmf)d八.防沖e 丑 Ee d Ebne和neb分別是場域單元數(shù)和邊界單元

5、數(shù)，E和Eb分別是單個場域單元積分區(qū)域和單個邊界單元積分區(qū)域。將基函數(shù)和權(quán)函數(shù)更換為單元形狀函數(shù)，在單元上將基函數(shù)和權(quán)函數(shù)用雙層下標(biāo)標(biāo)記：nenebKm,n= “aMmei ?Mne,j)dC- j (bM叫,Mn。,)d(17)e4 Ee=1 Eb電磁場數(shù)值計算結(jié)課大作業(yè)nenebKm,n= J(aVMmej NMne,j)dC J (bM%Mne)d(18)e 4 EejEb給定m和n的值，根據(jù)單元節(jié)點局部編號與整體編號的對應(yīng)關(guān)系，尋找單元上的i和j。單元形狀函數(shù)如下所示：nenebKm,n 八 、Ni J Nj2-q (bNjNj)d：(19)e=1 EeEbnenebRm C .(Ni

6、f)d八 (cNi)d】(20)e=Ee*得到得方程組以矩陣形式表示為：Ku = R(21)3、二維三角形有限元矩陣計算在二維三角形單元中(1)場域單元系數(shù)矩陣：Ke,i,j = .(a Ni .Nj)djE1/ 1 lyj書 _yj42aL(yi由- y42) (x也一x4OIxj 42 xj 書(2 )2(22)a ，(yi1-yi2)(yj1-yj2) ( -Xgj.2，其中，A為三角形單元面積。單元節(jié)點編號下標(biāo)采用1到3循環(huán)計數(shù)。(19)式中當(dāng)計算出的下標(biāo)大于 3時，取計算值減3為實際下標(biāo)。(2)邊界單元系數(shù)矩陣為1 :，22 1.Keb,i,j=-bNiNj d： e=b . X2

7、-小NiNjd(23)-e2令 l = (X2 -x1)2 (y2 -y1)2 ，得13 16Keb = -bl I ：:eb1/61/3(24)(3)場域單元右端項Q = E Nif d11(25)若源在單元上為常數(shù)，根據(jù)三角形單元解析積分公式:ARe:/一111心(26)否則根據(jù)源在單元節(jié)點上的函數(shù)值進行插值，這樣場域右端項為:電磁場數(shù)值計算結(jié)課大作業(yè)12 12 11-A12-eR112rj(27)(4)邊界右端項，可以表示為若c在單元上為常數(shù)，則, %(28)(29)若在單元上，c為線性插值函數(shù)，則(30)二、二維三角形單元有限元程序編寫及驗證1、程序思路開始采用PD具箱進行網(wǎng)格剖分，并

8、輸出剖分單元及節(jié)點等信息圖1二維三角形單元有限元法流程圖2、程序驗證電磁場數(shù)值計算結(jié)課大作業(yè)(1)算例1采用PDE工具箱給出的一個例子，對程序進行驗證。兩個圓形的金屬導(dǎo)體被放置在一張浸過鹽水的吸墨紙上， 吸墨紙起著平面薄導(dǎo)體的作用。這個問題的物理模型由拉普拉斯方程組成。求解場域以及邊界條件如圖2模型示意圖圖 2與式(31)所示。Z2u=0f(31)x = -1.2,-0.6 y0.6 。四=0, y = 0.6,-1.2x1.2 cnx = 1.2,-0.6y 0.6y = -0.6,-1.2x1.222u=1,(x+0.6)+ y =0.09u = _1,(x_0.6) +y2=0.09這里

9、令 。通過PDE工具箱進行網(wǎng)格音1J分，共得到512個單元，305個節(jié)點，如圖3所示。電磁場數(shù)值計算結(jié)課大作業(yè)-0.500.5圖3網(wǎng)格剖分結(jié)果使用PDE工具箱求解此問題，得到電壓等勢圖及電流場，如圖 4所示。Color: V Vector field: J-0.50.5圖4 PDE工具箱求得的電位及電流密度導(dǎo)出剖分所得單元、節(jié)點信息。使用自編程序進行求解，得到結(jié)果如下。電磁場數(shù)值計算結(jié)課大作業(yè)-0.50051圖5自編程序求得的電位及電流密度比較圖4 PDE工具箱求得的電位及電流密度圖4與圖5可知，自編程序得到了 PDE工具箱給出的結(jié)果。 圖上些許的差異，特別是等勢線與矢量場的不同，是由于自己繪

10、圖方式與PDE工具箱繪圖方式不同所導(dǎo)致。 進一步比較得到的電位值， 最大相對誤差為4.9982e-10, 可說明自編程序的正確性。(2)算例2在單位圓盤內(nèi)求解泊松方程：a，2u=P11(32)1 p =，+，22x yu = 0, x2 + y2 = 1網(wǎng)格剖分如圖6所示。電磁場數(shù)值計算結(jié)課大作業(yè)10.80.60.40.20-0,20.40.60.8-1-1-0.500.51圖6單位圓盤上的網(wǎng)格剖分通過PDE工具箱求解，得到電位分布如圖7所示。Numerical Solutionx101s8電磁場數(shù)值計算結(jié)課大作業(yè)導(dǎo)入網(wǎng)格后，用自編程序求解，得到結(jié)果如圖8所示。1F=f(node(1),nod

11、e(2),node(3);elseF=f*ones(1,3);endB=R(node(1),node(2),node(3).+delta/12*A*F.;R(node(1)=B(1);R(node(2)=B(2);R(node(3)=B(3);%R(node(1)=R(node(1)+delta*f/3;%R(node(2)=R(node(2)+delta*f/3;%R(node(3)=R(node(3)+delta*f/3;end%邊界單元for i=1:length(e(1,:)node(1)=e(1,i);node(2)=e(2,i);boundary=e(5,i);%所在的邊界編號if

12、 bc(1,boundary)=2%第二、第三類邊界條件b=-bc(2,boundary);%將pdeTool的公式轉(zhuǎn)化為書上的參數(shù)c=bc(3,boundary);%單元系數(shù)矩陣x(1)=p(1,node(1);y(1)=p(2,node(1);x(2)=p(1,node(2);y(2)=p(2,node(2);%長度l=sqrt( (x(1)-x(2)A2+(y(2)-y(1)A2);K(node(1),node(1)=K(node(1),node(1)+b*l/3;K(node(1),node(2)=K(node(1),node(2)+b*l/6;K(node(2),node(1)=K(

13、node(1),node(2)+b*l/6;K(node(2),node(2)=K(node(2),node(2)+b*l/3;%右端向量RR(node(1)=R(node(1)+c*l/2;R(node(2)=R(node(2)+c*l/2;else%第一類邊界條件u0=bc(3,boundary)/bc(2,boundary);% 求得給定的位函數(shù)值R=R-K(:,node(1)*u0;% (2)R(node(1)=u0;% (3)K(node(1),:)=0;% (4)K(:,node(1)=0;% (5)K(node(1),node(1)=1;%(1)% 注意課本上115頁這一語句在語

14、句(2)的前面.%這樣會使K矩陣對角線元素為零，使其不滿秩，無法求得正確的結(jié)果%猜測邊界上任何一點都會作為邊界元的開頭，不過為保險起見，先加入末端點%經(jīng)過驗證，確實不用下面這些%R=R-K(:,node(2)*u0;12電磁場數(shù)值計算結(jié)課大作業(yè)%R(node(2)=u0;%K(node(2),:)=0;%K(:,node(2)=0;%K(node(2),node(2)=1;endend%求解u1=KR;%最大誤差error=max(abs(u-u1)./u);%繪圖figure(1)x=linspace(min(p(1,:),max(p(1,:),20);y=linspace(min(p(2,:),max(p(2,:),20);xx,yy=meshgrid(x,y);uxy=tri2grid(p,t,u1,x,y);dx, dy=gradient(uxy);%繪制彩色三角形單元X,Y=meshgrid(p(1,:),p(