有限元綜合實驗報告

有限元方法及其程序設計綜合實驗報告目錄第一部分一．上機題目二．作業(yè)要求第二部分一．求偏分方程的變分形式二．區(qū)域剖分三．構(gòu)造有限元子空間四．建立有限元方程五．邊界條件處理六．求解代數(shù)方程組并輸出剛度矩陣第三部分一.輸出各節(jié)點處的計算結(jié)果二.圖形顯示結(jié)果，并做誤差分析比較附錄MATLAB有限元設計程序第一部分一、上機題目求解下列方程（1）其中求解區(qū)域，邊界，，，，，n表示單位外法向量。二、作業(yè)要求按照有限元方法的以下計算流程完成作業(yè)：求出方程的變分形式；區(qū)域剖分；構(gòu)造有限元子空間；建立有限元方程；邊界條件處理；求解代數(shù)方程組。第二部分一．求偏分方程的變分形式圖1求解區(qū)域及邊界示意圖令檢驗、試探函數(shù)空間分別為：利用Green公式，求方程（1）的變分形式求，使得（2）其中，，,,變分形式的計算過程如下：將原方程左右兩邊同時乘以，在內(nèi)積分，利用Green公式可得：左邊=====（3）右邊=（4）即=+=+ （5）二．區(qū)域剖分采用三角形單元對進行剖分。對方向均N等分（這里取N=4），計。h的選取必須保證把四個邊界的交點如（0,0）和（0,1）（1,0）（1,1）作為節(jié)點。首先對節(jié)點的總體編碼，（如圖2）。然后在每個單元上對節(jié)點進行局部編碼，生成總體編碼和局部編碼的對照表，如表2所示。對于節(jié)點，從邊界條件可以確定函數(shù)值的節(jié)點為Ⅰ類節(jié)點，個數(shù)為NE（NE=9）。未知函數(shù)值的節(jié)點為Ⅱ類，將對Ⅱ類每一個節(jié)點變量建立一個方程，故Ⅱ類節(jié)點的個數(shù)NF（即總節(jié)點數(shù)，本例中NF=16）是有限元空間的維數(shù)。對于Ⅱ類節(jié)點，確定影響元素集和影響節(jié)點集，如表1所示。圖2求解區(qū)域網(wǎng)格剖分示意圖表1類節(jié)點、影響元素集、影響節(jié)點集對照表（節(jié)點）（影響元素集）（影響節(jié)點集）72349101123678111284561112133478912139678131415458910131410815165910141512101112171819781112131617131213141920218912131417181414151621222391013141518191516232410141519201718192025262712131617182122182021222728291314171819222319222324293031141518192023242024313215192024252226272817182122232328293018192223242430313219202324252532202425表2點局部編碼與整體編碼對照表一二三四五六七八九十十一十二十三十四十五十六Ⅰ172839410612713814915Ⅱ263748597118129131014Ⅲ627384951171281391410十七十八十九二十二一二二二三二四二五二六二七二八二九三十三一三二Ⅰ11171218131914201622172318241925Ⅱ12161317141815191721182219232024Ⅲ16121713181419152117221823192420用MATLAB程序?qū)崿F(xiàn)局部編碼和整體編碼對照表的方法ⅡⅢ從表2中可以發(fā)現(xiàn)，在32(2×N^2)個單元中，局部編碼為Ⅰ的節(jié)點，其整體編碼可以按照單元數(shù)為奇數(shù)和偶數(shù)分成兩類，即奇數(shù)單元其總體編碼為1、2、3、4；6、7、8、9；11、12、13、14；16、17、18、19；偶數(shù)單元其總體編碼為2、3、4、5；7、8、9、10；12、13、14、15；17、18、19、20。這里的分號表示把其中的4（N）個數(shù)放在一組，可以組成4×4（N×N）的矩陣?？梢园l(fā)現(xiàn)，如果生成一個5×5（N+1×N+1）矩陣，其中的元素從1到25（N+1^2），先按行再按列的原則分布。很容易發(fā)現(xiàn)，奇數(shù)單元局部編碼為Ⅰ的節(jié)點，其整體編碼矩陣為5×5（N+1×N+1）矩陣左上角4×4（N×N）子矩陣；同理，偶數(shù)單元局部編碼為Ⅰ的節(jié)點，其整體編碼矩陣為5×5（N+1×N+1）矩陣右下角4×4（N×N）子矩陣。同時在表2中可以發(fā)現(xiàn)，奇數(shù)單元局部編碼為Ⅰ的節(jié)點與局部編碼為Ⅱ、Ⅲ的節(jié)點的總體編碼的差值分別為1和5，由此可以根據(jù)上面求的局部編碼為Ⅰ節(jié)點的總體編碼值計算局部編碼為Ⅱ、Ⅲ的節(jié)點的總體編碼；同理，偶數(shù)單元局部編碼為Ⅰ的節(jié)點與局部編碼為Ⅱ、Ⅲ的節(jié)點的總體編碼的差值分別為-1和-5，由此，可以根據(jù)上面求的局部編碼為Ⅰ節(jié)點的總體編碼值計算局部編碼為Ⅱ、Ⅲ的節(jié)點的總體編碼。至此，所有單元的局部編碼和整體編碼對照關(guān)系可以確定下來，以下為MATLAB程序的實現(xiàn)三．構(gòu)造有限元子空間建立有限元子空間，其中=，NF=16。節(jié)點基是由分片多項式拼成的。且和，則（2）式離散形式為求，使得，（6）其中，，。選取三角形單元的自由度集合是：由可知：，表示的值。四．建立有限元方程由式（6）可知：，其中，，，可以得到=+=+因為上式對任意均成立，所以上式可以簡化為=+，=+，+，所以可得：=+，（7）需要對NF個節(jié)點（這里NF=16）分別建立代數(shù)方程。內(nèi)部節(jié)點代數(shù)方程建立，以第7個方程為例（j=7）圖3節(jié)點7的影響元素集和影響節(jié)點集示意圖如圖3所示，找出節(jié)點7的影響元素集和影響節(jié)點集，即把j=7代入式（7）中得到（8）節(jié)點基是由，，給出的。=2==其中是標準元。是由，，給出的，即類似地=2+2+2=同理，可以求得其他系數(shù)其中是節(jié)點的影響元素集的元素個數(shù)。綜上所述，第7個方程為：=邊界上節(jié)點代數(shù)方程建立，以第23個方程為例（j=23）邊界節(jié)點代數(shù)方程建立原則和規(guī)定：由待求方程和圖4可以得出，在邊界上節(jié)點5、10、15、20、25的線積分有值，對節(jié)點線積分的積分長度采取左移的原則，如圖4中的箭頭所示方向（例如對節(jié)點23求線積分時，積分長度為節(jié)點22和23之間的長度區(qū)間）；在邊界上節(jié)點22、23、24、25的線積分有值，對節(jié)點線積分的積分長度采取下移的原則，如圖4中的箭頭所示方向，（例如對節(jié)點15求線積分時，積分長度為節(jié)點15和10之間的長度區(qū)間）。其他形式均與內(nèi)部節(jié)點相同。圖4邊界節(jié)點線積分積分區(qū)間示意圖圖5節(jié)點23的影響元素集合影響節(jié)點集示意圖如圖5所示，找出節(jié)點23的影響元素集和影響節(jié)點集，即把j=23代入式（7）中得到（9）類似于上述方程求解過程，可以求得各系數(shù)為：，，其中在第28個單元中綜上所述，第23個方程是邊界上節(jié)點代數(shù)方程建立，以第15個方程為例（j=15）圖6節(jié)點23的影響元素集合影響節(jié)點集示意圖如圖6所示，找出節(jié)點15的影響元素集和影響節(jié)點集，即把j=15代入式（7）中得到（10）類似于上述方程求解過程，可以求得各系數(shù)為：，，其中在第16個單元中綜上所述，第15個方程是五．邊界條件處理由于在上述MATLAB實現(xiàn)的過程中，直接把25（N+1^2）個節(jié)點的都考慮進去了，但是實際求取中，要刪除第Ⅰ類（已知值）節(jié)點所在剛度矩陣的行和列（例如本例中，對于節(jié)點6，就需要把剛度矩陣的第六行和第六列刪除）。同理，右端列陣也相應要刪除已知節(jié)點所在的行，對于第Ⅱ類節(jié)點的方程中含有第Ⅰ類節(jié)點的運算，計算相應結(jié)果，把它放在方程的右端。本例中已知節(jié)點的值為0，故第Ⅱ類節(jié)點的方程中不含有第Ⅰ類節(jié)點的運算，給計算帶了了方便。經(jīng)過邊界條件處理后得到的方程為：六．求解代數(shù)方程組并輸出剛度矩陣由此得到NF（這里NF=16）個未知量的NF各方程組則有表3剛度矩陣系數(shù)K（N=4）第三部分一.輸出各個節(jié)點處的計算結(jié)果1．N=4時各節(jié)點的函數(shù)值如表4所示表4N=4時的計算結(jié)果坐標00.250.50.751.00000000.2500.05500.11370.18130.27060.500.10710.21330.32280.43810.7500.16440.31020.44790.57660.100.24580.41660.56790.7353二.圖形顯示結(jié)果，并做誤差分析比較圖7N=4時的計算結(jié)果圖7N=8時的計算結(jié)果圖8N=16時的計算結(jié)果表5不同密度網(wǎng)格下計算結(jié)果的對比每段邊劃分的段數(shù)（N）1,1）點處的計算函數(shù)值相鄰N值之間的相對誤差40.735380.67098.76﹪120.64623.68﹪160.63282.07﹪附錄：%---------------------------------------------clear%-----------------已知參量說明----------------N=input('PleaseinputN:\n=');h=1/N;%-----------生成N等分局部編碼和整體編碼對照表---------fori=1:(N+1)forj=1:(N+1)a(i,j)=(i-1)*(N+1)+j;endendm=a(1:N,1:N);%提取奇數(shù)單元第一類節(jié)點總體編碼矩陣n=a(2:N+1,2:N+1);%提取偶數(shù)單元第一類節(jié)點總體編碼矩陣II1=[];II2=[];fori=1:NII1=[II1m(i,:)];%合成奇數(shù)單元第一類節(jié)點總體編碼行向量II2=[II2n(i,:)];%合成偶數(shù)單元第一類節(jié)點總體編碼行向量endforLEE=1:2*N^2%合并奇、偶單元第一類節(jié)點總體編碼行向量ifmod(LEE,2)==1II(1,LEE)=II1((LEE+1)/2);elseII(1,LEE)=II2(LEE/2);endendforLEE=1:2:2*N^2-1%計算奇數(shù)單元第一類、第二類、第三類節(jié)點總體編碼II(2,LEE)=II(1,LEE)+1;II(3,LEE)=II(1,LEE)+(N+1);endforLEE=2:2:2*N^2%計算偶數(shù)單元第一類、第二類、第三類節(jié)點總體編碼II(2,LEE)=II(1,LEE)-1;II(3,LEE)=II(1,LEE)-(N+1);end%------------產(chǎn)生對稱剛度舉矩陣程序---------------------k=[1+h^2/12-1/2-1/2;-1/21/2+h^2/120;-1/201/2+h^2/12]+h^2/12;%單元剛度矩陣K=zeros((N+1)^2);%總剛度矩陣初始化fori=1:(N+1)^2forj=1:(N+1)^2forLEE=1:2*N^2P1=0;%記錄節(jié)點i的局部編碼P2=0;%記錄節(jié)點j的局部編碼forND=1:3ifII(ND,LEE)==iP1=ND;breakendendforND=1:3ifII(ND,LEE)==jP2=ND;breakendendif(P1&P2)==1K(i,j)=K(i,j)+k(P1,P2);endendendend%-----------------------------求(f,v)-----------------------------L(1:(N+1)^2)=0;fori=1:(N+1)^2forLEE=1:2*N^2forND=1:3ifII(ND,LEE)==iL(i)=L(i)+1;%計算節(jié)點i所在的單元數(shù)endendendendfori=1:(N+1)^2ifmod(i,(N+1))~=0x=(mod(i,(N+1))-1)*h;y=fix(i/(N+1))*h;f(i)=(x^2*cos(y)+y^2*sin(x)+x^2*y^2)*L(i)*h^2/6;elsex=1;y=(i/(N+1)-1)*h;f(i)=(x^2*cos(y)+y^2*sin(x)+x^2*y^2)*L(i)*h^2/6;endend%------------求邊界上2上的函數(shù)值g2---------------------symsy;fori=2*(N+1):(N+1):(N+1)^2yy=(i/(N+1)-1)*h;a=[1y1;1-hyy1;1yy-h1];p(i)=(1/h^2)*det(a);g2(i)=int(p(i)*cos(y),yy-h,yy);%求邊界2上的線積分end%-------------求邊界上3上的函數(shù)值g3--------------------symsx;fori=(N+1)^2-N+1:(N+1)^2ifmod(i,N+1)~=0xx=(mod(i,N+1)-1)*h;a=[x11;xx-h11;xx1-h1];p(i)=(1/h^2)*det(a);g3(i)=int(p(i)*cos(x),xx-h,xx);%求邊界3上的線積分elsexx=1;a=[x11;xx-h11;xx1-h1];p(i)=(1/h^2)*det(a);g3(i)=int(p(i)*cos(x),xx-h,xx);%求邊界3上的線積分endend%------------