有限元的matlab编程

输出求解结果：
输出位移
fprintf( '节点位移\n' ) ;
for i=1:1:10 disp(['节点号',num2str(i),' x方向位移：',num2str(u(2*i-1,1)),' y方 向位移：',num2str(u(2*i,1))]); end
输出节点力
fprintf( '\n\n节点力\n' ) ; for ie=1:1:21 nodef=NodeForce( ie ); disp( ['单元号: ',num2str(ie),' 节点号：',num2str(Element(ie,1)),'
用户自定义网架（网架信息的录入，包括节点、单元、截面、弹性模量等）
if e==0%选择自定义网架
Node=input(‘定义节点编号及对应坐标，按[1 x1 y1 z1;2 x2 y2 z2;...]输入’);%形成节点储存矩阵
Men=input(‘定义单元与节点的关系，按[1 node1 node2;2 node3 node4;...]输入,node1<node2, 依次类推’);%形成单元储存矩阵 Msum=length(Men);%查找网架录入的单元数
正放四角锥网架定义
if e==1
定义网架上层节点
hu=input('输入网架上层节点行数'); %定义网架上层节点的行数
lu=input('输入网架上层节点列数'); %定义网架上层节点的列数
dis_xu=input('输入网架上层节点列间距'); %定义网架上层的行间距
dis_yu=input('输入网架上层节点行间距'); %定义网架上层的列间距
网架上层横向单元的拓扑
for i=1:lu
for j=1:hu-1
Men((i-1)*(hu-1)+j+(lu-1)*hu,2)=(j-1)*lu+i; Men((i-1)*(hu-1)+j+(lu-1)*hu,3)=j*lu+i; end end
Men((i-1)*(lu-1)+j,3)=(i-1)*lu+j+1;
节点号：',num2str(Element(ie,2)),' 轴力:',num2str(nodef(1))] ) ;
end
例二：网架
思路分析
几点说明
网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成 网架的基本单元有三角锥，三棱体，正方体，截头四角锥等。鉴于网架的形
式较多，本程序提供一种通用的网架输入方法，但录入较为繁琐，同时提供
0
E*A/L 0
0
0 0];
T = TransformMatrix( ie ) ;
k = T*k*transpose(T) ;% transpose(T) 为T的转置矩阵2
集成整体刚度矩阵K
K=zeros(20,20);%用来存储整体刚度矩阵 在下面的集成中，将总刚看成10*10的矩阵，每个元素为2*2的小矩阵
集成总刚的非对角线元素(这里的元素指2*2的小矩阵）
for ie=1:1:21 %按单元顺序进行循环
k=PlaneTrussElementStiffness(ie); %计算第ie个单元的单刚 m=Element(ie,1); %ie单元的首节点号
n=Element(ie,2); %ie单元的末节点号
K(2*m,2*n-1)=k(2,3);
K(2*m,2*n)=k(2,4);
集成总刚的对角线元素（这里的元素指2*2的小矩阵）
for i=1:1:10 %按节点的顺序循环 for j=1:1:21 %对于每个节点，再按单元的顺序循环 k=PlaneTrussElementStiffness(j); if Element(j,1)==I %如果i节点为j单元的首节点 K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(1,1); K(2*i-1,2*i)=K(2*i-1,2*i)+k(1,2); K(2*i,2*i-1)=K(2*i,2*i-1)+k(2,1); K(2*i,2*i)=K(2*i,2*i)+k(2,2); end if Element(j,2)==i %如果i节点为j单元的末节点 K(2*i-1,2*i-1)=K(2*i-1,2*i-1)+k(3,3); K(2*i-1,2*i)=K(2*i-1,2*i)+k(3,4); K(2*i,2*i-1)=K(2*i,2*i-1)+k(4,3); K(2*i,2*i)=K(2*i,2*i)+k(4,4); end end end
求解
f=f*1e15;
u=K\f;
求解轴力：
获取单元两端的节点号
i = Element( ie, 1 ) ;%ie为单元号 j = Element( ie, 2 ) ;
获取单元两端的节点位移
uElement = zeros( 4, 1 ) ; uElement( 1:2 ) = u( (i-1)*2+1:(i-1)*2+2 ) ; uElement( 3:4 ) = u( (j-1)*2+1:(j-1)*2+2 ) ;
求解位移：
u=zeros(20);
根据约束情况修改总刚，采用对角元素置1法
for i=1:1:20 K(1,i)=0; K(2,i)=0; K(18,i)=0; K(i,1)=0;K(i,2)=0; K(i,18)=0; end %自由度1、2、18被约束了,所在的行和列的其他元素都改为0 K(1,1)=1;%对角线元素置1 K(2,2)=1; K(18,18)=1; K=K*1e15;%乘以一个大数，减小计算误差
Element=zeros(21,2); for i=1:2:7 Element(5/2*i-3/2,:)=[i,i+1]; Element(5/2*i-1/2,:)=[i,i+2]; Element(5/2*i+1/2,:)=[i,i+3]; end
for i=2:2:8
Element(5*i/2-1,:)=[i,i+1]; Element(5*i/2,:)=[i,i+2];
Node((i-1)*ld+j+hu*lu,4)=0;
end end
节点编号的录入
Nsum=length(Node); %查询网架的节点数 for i=1:Nsum %将节点编号录入节点矩阵 Node(i,1)=i; end
网架上层纵向单元的拓扑
for i=1:hu
for j=1:lu-1 Men((i-1)*(lu-1)+j,2)=(i-1)*lu+j;
K(2*m-1,2*n-1)=k(1,3); K(2*m-1,2*n)=k(1,4);
m=Element(ie,2); %ie单元的末节点号 n=Element(ie,1); %ie单元的首节点号 K(2*m-1,2*n-1)=k(3,1); K(2*m-1,2*n)=k(3,2); K(2*m,2*n-1)=k(4,1); K(2*m,2*n)=k(4,2); end
有限元编程示例
例一：桁架
题目描述:
如下图所示的平面桁架，杆件长度、弹性模量、截面积以 及所受节点力P的大小可以自行定义。求节点位移及杆件轴 力。
解题思路：
• 建立模型 • 集成总刚 • 求解位移
• 求解杆件轴力
• 输出结果
建立模型：
模型相关参数输入
H=input('竖杆长度(m):'); L=input('水平杆长度(m):'); E=input('杆件弹性模量(Gpa):'); A=input('杆件截面积(m^2):'); a=input('节点力P(kN)：');
function T = TransformMatrix( ie )%ie为单元号
c = (xj-xi)/L ;
s = (yj-yi)/L ;
T=[ c -s s c 0 0 0 0
0 0
0 0
c -s
s c ];
计算单元刚度矩阵k
k = [ E*A/L 0 -E*A/L 0
0
-E*A/L
0
0 0
计算单元的节点力
k = PlaneTrussElementStiffness( ie ) ; nodef = k *uElement ;%整体坐标下的节点力
T = TransformMatrix( ie ) ;%计算坐标转换矩阵
nodef = transpose( T ) * nodef ;%从整体坐标转换到局部坐标
yi = Node( Element( ie, 1 ), 2 ) ;
xj = Node( Element( ie, 2 ), 1 ) ; yj = Node( Element( ie, 2 ), 2 ) ;
计算杆件长度
L = ( (xj-xi)^2 + (yj-yi)^2 )^(1/2) ;
计算从局部坐标到整体坐标的坐标转换矩阵T
定义单元属性的输入方式
Cont1=input('定义单元实常数，若所有杆件截面面积和弹性模量不变，则输入0，否则输入1');
单元属性相同 if Cont1==0 AE1=input('请输入统一的截面面积与弹性模量，按[A E]输入'); AE=zeros(Msum,3); AE(:,1)=1:Msum;AE(:,2)=AE1(1,1);AE(:,3)=AE1(1,2); 单元属性不同 else AE=input('请输入相应单元的截面面积与弹性模量，按[1,A1 E1;2,A2 E2;...]输入'); end