有限元方法及其程序设计
综合实验报告
目录
第一部分
一．上机题目
二．作业要求
第二部分
一．求偏分方程的变分形式
二．区域剖分
三．构造有限元子空间
四．建立有限元方程
五．边界条件处理
六．求解代数方程组并输出刚度矩阵
第三部分
一.输出各节点处的计算结果
二.图形显示结果，并做误差分析比较
附录
MATLAB有限元设计程序
第一部分
一、上机题目
求解下列方程
其中求解区域
，边界
，
，
，
，
，n表示单位外法向量。
二、作业要求
按照有限元方法的以下计算流程完成作业：
1. 求出方程的变分形式；
2. 区域剖分；
3. 构造有限元子空间；
4. 建立有限元方程；
5. 边界条件处理；
6. 求解代数方程组。
第二部分
一．求偏分方程的变分形式
图1 求解区域及边界示意图
令检验、试探函数空间分别为：
利用Green公式，求方程（1）的变分形式
求
，使得
（2）
其中，
，
, 
, 
变分形式的计算过程如下：
将原方程左右两边同时乘以
，在
内积分，利用Green公式可得：
左边=
=
=
=
=
（3）
右边=
（4）
即 
=+
=+
（5）
二．区域剖分
采用三角形单元对
进行剖分。对
方向均N等分（这里取N=4），计
。h的选取必须保证把
四个边界的交点如（0,0）和（0,1）（1,0）（1,1）作为节点。首先对节点
的总体编码
，
（如图2）。然后在每个单元上对节点进行局部编码，生成总体编码和局部编码的对照表，如表2所示。
对于节点，从边界条件可以确定函数值的节点为Ⅰ类节点，个数为NE（NE=9）。未知函数值的节点为Ⅱ类，将对Ⅱ类每一个节点变量建立一个方程，故Ⅱ类节点的个数NF（即总节点数
，本例中NF=16）是有限元空间的维数。对于Ⅱ类节点，确定影响元素集和影响节点集，
如表1所示。
图2 求解区域网格剖分示意图
表1 类节点、影响元素集、影响节点集对照表
表2 点局部编码与整体编码对照表
用MATLAB程序实现局部编码和整体编码对照表的方法ⅡⅢ
从表2中可以发现，在32(2×N^2)个单元中，局部编码为Ⅰ的节点，其整体编码可以按照单元数为奇数和偶数分成两类，即奇数单元其总体编码为1、2、3、4；6、7、8、9；11、12、13、14；16、17、18、19；偶数单元其总体编码为2、3、4、5；7、8、9、10；12、13、14、15；17、18、19、20。这里的分号表示把其中的4（N）个数放在一组，可以组成4×4（N×N）的矩阵。可以发现，如果生成一个5×5（N+1×N+1）矩阵，其中的元素从1到25（N+1^2），先按行再按列的原则分布。很容易发现，奇数单元局部编码为Ⅰ的节点，其整体编码矩阵为5×5（N+1×N+1）矩阵左上角4×4（N×N）子矩阵；同理，偶数单元局部编码为Ⅰ的节点，其整体编码矩阵为5×5（N+1×N+1）矩阵右下角4×4（N×N）子矩阵。
同时在表2中可以发现，奇数单元局部编码为Ⅰ的节点与局部编码为Ⅱ、Ⅲ的节点的总体编码的差值分别为1和5，由此可以根据上面求的局部编码为Ⅰ节点的总体编码值计算局部编码为Ⅱ、Ⅲ的节点的总体编码；同理，偶数单元局部编码为Ⅰ的节点与局部编码为Ⅱ、Ⅲ的节点的总体编码的差值分别为-1和-5，由此，可以根据上面求的局部编码为Ⅰ节点的总体编码值计算局部编码为Ⅱ、Ⅲ的节点的总体编码。至此，所有单元的局部编码和整体编码对照关系可以确定下来，以下为MATLAB程序的实现
三．构造有限元子空间
建立有限元子空间
，其中
=
，NF=16。节点基
是由分片多项式拼成的。且
和
，则（2）式离散形式为
求
，使得
， 
（6）
其中，
，
。
选取三角形单元的自由度集合是：
由可知：
，
表示
的值。
四．建立有限元方程
由式（6）可知：
，
其中，
，，可以得到
=
+
=
+
因为上式对任意
均成立，所以上式可以简化为
=
+
， 
=+，
+，
所以可得：
=
+， （7）
需要对NF个节点（这里NF=16）分别建立代数方程。
1）内部节点代数方程建立，以第7个方程为例（j=7）
图3 节点7的影响元素集和影响节点集示意图
如图3所示，找出节点7的影响元素集和影响节点集，即把j=7代入式（7）中得到
（8）
节点基是由
，
，
给出的。
=2
=
=
其中
是标准元。是由，，给出的，即
类似地
=2
+2
+2
=
同理，可以求得其他系数
其中
是节点
的影响元素集
的元素个数。
综上所述，第7个方程为：
=
2） 边界
上节点代数方程建立，以第23个方程为例（j=23）
边界节点代数方程建立原则和规定：
由待求方程和图4可以得出，在边界
上节点5、10、15、20、25的线积分
有值，对节点线积分的积分长度采取左移的原则，如图4中的箭头所示方向（例如对节点23求线积分时，积分长度为节点22和23之间的长度区间）；在边界上节点22、23、24、25的线积分
有值，对节点线积分的积分长度采取下移的原则，如图4中的箭头所示方向，（例如对节点15求线积分时，积分长度为节点15和10之间的长度区间）。其他形式均与内部节点相同。
图4 边界节点线积分积分区间示意图
图5 节点23的影响元素集合影响节点集示意图
如图5所示，找出节点23的影响元素集和影响节点集，即把j=23代入式（7）中得到
（9）
类似于上述方程求解过程，可以求得各系数为：
，
，
其中
在第28个单元中
综上所述，第23个方程是
3） 边界上节点代数方程建立，以第15个方程为例（j=15）
图6 节点23的影响元素集合影响节点集示意图
如图6所示，找出节点15的影响元素集和影响节点集，即把j=15代入式（7）中得到
（10）
类似于上述方程求解过程，可以求得各系数为：
，
，
其中
在第16个单元中
综上所述，第15个方程是
五．边界条件处理
由于在上述MATLAB实现的过程中，直接把25（N+1^2）个节点的都考虑进去了，但是实际求取中，要删除第Ⅰ类（已知值）节点所在刚度矩阵的行和列（例如本例中，对于节点6，就需要把刚度矩阵的第六行和第六列删除）。同理，右端列阵也相应要删除已知节点所在的行，对于第Ⅱ类节点的方程中含有第Ⅰ类节点的运算，计算相应结果，把它放在方程的右端。
本例中已知节点的值为0，故第Ⅱ类节点的方程中不含有第Ⅰ类节点的运算，给计算带了了方便。
经过边界条件处理后得到的方程为：
六．求解代数方程组并输出刚度矩阵
由此得到NF（这里NF=16）个未知量
的NF各方程组
则有
表3 刚度矩阵系数K（N=4）
第三部分
一.输出各个节点处的计算结果
1．N=4时各节点的函数值如表4所示
表4 N=4时的计算结果
二.图形显示结果，并做误差分析比较
图7 N=4时的计算结果
图7 N=8时的计算结果
图8 N=16时的计算结果
表5 不同密度网格下计算结果的对比
附录：
%---------------------------------------------
clear
%-----------------已知参量说明----------------
N=input('Please in put N:\n =');
h=1/N;
%-----------生成N等分局部编码和整体编码对照表---------
for i=1:(N+1)
for j=1:(N+1)
a(i,j)=(i-1)*(N+1)+j;
end
end
m=a(1:N,1:N);%提取奇数单元第一类节点总体编码矩阵
n=a(2:N+1,2:N+1);%提取偶数单元第一类节点总体编码矩阵
II1=[];
II2=[];
for i=1:N
II1=[II1 m(i,:)];%合成奇数单元第一类节点总体编码行向量
II2=[II2 n(i,:)];%合成偶数单元第一类节点总体编码行向量
end
for LEE=1:2*N^2 %合并奇、偶单元第一类节点总体编码行向量
if mod(LEE,2)==1
II(1,LEE)=II1((LEE+1)/2);
else
II(1,LEE)=II2(LEE/2);
end
end
for LEE=1:2:2*N^2-1 %计算奇数单元第一类、第二类、第三类节点总体编码
II(2,LEE)=II(1,LEE)+1;
II(3,LEE)=II(1,LEE)+(N+1);
end
for LEE=2:2:2*N^2 %计算偶数单元第一类、第二类、第三类节点总体编码
II(2,LEE)=II(1,LEE)-1;
II(3,LEE)=II(1,LEE)-(N+1);
end
%------------产生对称刚度举矩阵程序---------------------
k=[1+h^2/12 -1/2 -1/2;-1/2 1/2+h^2/12 0;-1/2 0 1/2+h^2/12]+h^2/12;%单元刚度矩阵
K=zeros((N+1)^2);%总刚度矩阵初始化
for i=1:(N+1)^2
for j=1:(N+1)^2
for LEE=1:2*N^2
P1=0;%记录节点i的局部编码
P2=0;%记录节点j的局部编码
for ND=1:3
if II(ND,LEE)==i
P1=ND;
break
end
end
for ND=1:3
if II(ND,LEE)==j
P2=ND;
break
end
end
if (P1&P2)==1
K(i,j)=K(i,j)+k(P1,P2);
end
end
end
end
%-----------------------------求(f,v)-----------------------------
L(1:(N+1)^2)=0;
for i=1:(N+1)^2
for LEE=1:2*N^2
for ND=1:3
if II(ND,LEE)==i
L(i)=L(i)+1; %计算节点i所在的单元数
end
end
end
end
for i=1:(N+1)^2
if mod(i,(N+1))~=0
x=(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;
else
x=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;
end
end
%------------求边界上2上的函数值g2---------------------
syms y;
for i=2*(N+1):(N+1):(N+1)^2
yy=(i/(N+1)-1)*h;
a=[1 y 1;1-h yy 1;1 yy-h 1];
p(i)=(1/h^2)*det(a);
g2(i)=int(p(i)*cos(y),yy-h,yy);%求边界2上的线积分
end
%-------------求边界上3上的函数值g3--------------------
syms x;
for i=(N+1)^2-N+1:(N+1)^2
if mod(i,N+1)~=0
xx=(mod(i,N+1)-1)*h;
a=[x 1 1;xx-h 1 1;xx 1-h 1];
p(i)=(1/h^2)*det(a);
g3(i)=int(p(i)*cos(x),xx-h,xx);%求边界3上的线积分
else
xx=1;
a=[x 1 1;xx-h 1 1;xx 1-h 1];
p(i)=(1/h^2)*det(a);
g3(i)=int(p(i)*cos(x),xx-h,xx);%求边界3上的线积分
end
end
%----------------求右端列阵F的合成----------------------
for i=1:(N+1)^2
F(i)=f(i)+g2(i)+g3(i);
end
%---------------提取第二类节点刚度矩阵------------------
K(1:(N+1),:)=[];%删除边界1上节点所在刚度矩阵的所有行
K(:,1:(N+1))=[];%删除边界1上节点所在刚度矩阵的所有列
for i=1:N
K((i-1)*N+1,:)=[];%删除边界4上节点所在刚度矩阵的所有行
K(:,(i-1)*N+1)=[];%删除边界4上节点所在刚度矩阵的所有列
end
%----------------提取第二类节点右端列阵F-----------------
F=F';
F(1:(N+1),:)=[];%删除边界1上节点所在刚度矩阵的所有行
for i=1:N
F((i-1)*N+1,:)=[];%删除边界4上节点所在刚度矩阵的所有行
end
%---------------求解代数方程组---------------------
U=inv(K)*F;
U=eval(U);%让输出结果由分数转化为小数
%---------------图形显示结果---------------------
U=U';
Z(1:(N+1))=0;
for i=1:N
Z=[Z;0 U(1,N*i-N+1:N*i)];%添加第一类节点的值，使Z矩阵能表达所有节点的值
end
x=0:h:1;
y=0:h:1;
axis([0 1 0 1]);
surf(x,y,Z)
%---------------------------------------------------