3结点三角形单元有限元程序MATLAB语言

有限元课作业,我还是MATLAB稍微会点,C好久不用就没用C写,这个是大作业报告。

3结点三角形单元有限元程序（MATLAB语言）

该程序包括以下6个部分：

1．主程序tri_fem：用于数据的录入和其他程序的调用； 2．总刚程序Kf：计算结构的总体刚度；

3．各结点位移求解程序xf：求解各结点的位移；

4．线性方程组求解程序Jordan：Gauss-Jordan法求解非约束结点的位移；

5．应力应变程序ss：由各结点位移求解各单元内的三个结点的应力stress和应变strain；

6．数据录入程序input：录入材料、几何尺寸、单元编号和结点编号、位移约束和已知载荷等。

以课本P25页例2.2为例，其input程序为

function [E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input()

E=2.1e11; v=1/3; t=1; %杨氏模量Pa，泊松比，厚度 EN=2; %单元数

ecode=[3 1 2; %单元编号 单元1 3-1-2;单元2 1-3-4 1 3 4];

NN=4; %结点数 node=[0 0; %各结点坐标 2 0; 2 1;

0 1];

RN=2; %被约束的位移数 RC=[1 4]; %被约束的结点 PN=2; %有载荷的结点数

%PC（1）表示有载荷的结点，PC(2)表示各结点的力，PC（3）表示载荷方向，0为x方向， 1为y方向

PC=[2 3; -1/2 -1/2;

1 1]; %结点2、3分别有y负方向上-1/2N的力作用

有限元课作业,我还是MATLAB稍微会点,C好久不用就没用C写,这个是大作业报告。

在matlab环境下，输入则程序运行结果如下：

该程序求解的结点位移结果和结点应力结果与课本给出的结果一致。

有限元课作业,我还是MATLAB稍微会点,C好久不用就没用C写,这个是大作业报告。

附录：

%%-------------平面三角形单元有限元法---------------------%% function [x strain stress]=tri_fem()

[E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input; %调入已定材料、几何尺寸以及单元和结点编号及约束和载荷分布 [n m]=size(ecode); if EN~=n

error('Wrong elementnumber EN or wrong elementcode ecode!'); return; end

[n m]=size(node); if NN~=n

error('Wrong nodenumber NN or wrong node-coordinate node!'); return; end

e=zeros(EN,6);

A=zeros(EN,1); %面积 for i=1:EN

e(i,:)=[node(ecode(i,1),:),node(ecode(i,2),:),node(ecode(i,3),:)]; %各三角形单元的节点坐标

D=[1,node(ecode(i,1),:) 1,node(ecode(i,2),:) 1,node(ecode(i,3),:)]; A(i,1)=1/2*det(D); end

%% 形成单元参数 b=zeros(EN,3);

c=zeros(EN,3); %各单元参数初始化 for i=1:EN

b(i,1)=e(i,4)-e(i,6); b(i,2)=e(i,6)-e(i,2); b(i,3)=e(i,2)-e(i,4); c(i,1)=-e(i,3)+e(i,5); c(i,2)=-e(i,5)+e(i,1); c(i,3)=-e(i,1)+e(i,3); end

%% 求得总刚，并引入约束和载荷求得各结点位移

K=Kf(E,v,t,EN,ecode,NN,A,b,c); %调用函数Kf，求得结构的总体刚度矩阵 x=xf(NN,RN,RC,PN,PC,K); %调用函数xf,求得各结点位移 %% 求解应力应变

[strain stress]=ss(E,v,EN,ecode,A,b,c,x);

%% 单元刚度矩阵与结构刚度矩阵 function K=Kf(E,v,t,EN,ecode,NN,A,b,c)

Ke=zeros(6,6); %单元的刚度矩阵，初始为6*6阶零矩阵

有限元课作业,我还是MATLAB稍微会点,C好久不用就没用C写,这个是大作业报告。

K=zeros(NN*2,NN*2); %结构的总体刚度矩阵，初始为零矩阵 for m=1:1:EN %m为单元号

for i=1:1:3 for j=1:1:3

Ke(2*i-1,2*j-1)=b(m,i)*b(m,j)+(1-v)*c(m,i)*c(m,j)/2; Ke(2*i-1,2*j)=v*c(m,i)*b(m,j)+(1-v)*b(m,i)*c(m,j)/2; Ke(2*i,2*j-1)=v*b(m,i)*c(m,j)+(1-v)*c(m,i)*b(m,j)/2; Ke(2*i,2*j)=c(m,i)*c(m,j)+(1-v)*b(m,i)*b(m,j)/2; end end

Ke=E*t/(4*(1-v^2)*A(EN)).*Ke; %获得单元m的刚度矩阵

Kb=mat2cell(Ke,ones(1,3)*2,ones(1,3)*2); %将单元矩阵Ke分为3*3块 set1=ones(1,NN)*2;

Ka=mat2cell(zeros(NN*2,NN*2),set1,set1); %将总刚K分为NN*NN块 for i=1:1:3 for j=1:1:3

Ka(ecode(m,i),ecode(m,j))=Kb(i,j); %各单元刚度矩阵整体编号,并叠加 end end

K=K+cell2mat(Ka); end

%分块矩阵K合成一个矩阵K

%% 引入位移向量和右端项

function x=xf(NN,RN,RC,PN,PC,K)

x=ones(NN*2,1); %位移初始为0向量 for i=1:RN

x(RC(i)*2-1)=0; x(RC(i)*2)=0;

end %被约束结点位移为0 %%----------------引入已知结点载荷-------------% px=zeros(NN*2,1); %载荷初始为0向量 for i=1:PN

if PC(3,i)==1

px(PC(1,i)*2)=PC(2,i); else if PC(3,i)==0

px(PC(1,i)*2-1)=PC(2,i); end end end

%%----------------引入已知结点载荷-------------%

%%----------------求解非约束结点的位移X-------------% set1=ones(1,NN)*2;

有限元课作业,我还是MATLAB稍微会点,C好久不用就没用C写,这个是大作业报告。

Ka=mat2cell(K,set1,set1); pxa=mat2cell(px,set1,1);

AN=zeros(2*(NN-RN),2*(NN-RN));

ANa=mat2cell(AN,ones(1,NN-RN)*2,ones(1,NN-RN)*2); bn=zeros(2*(NN-RN),1);

bna=mat2cell(bn,ones(1,NN-RN)*2,1); BN=zeros(2*RN,2*(NN-RN));

BNa=mat2cell(BN,ones(1,RN)*2,ones(1,NN-RN)*2); m=1;

for i=1:1:NN if x(2*i)==1 M(m)=i; m=m+1; end end

for i=1:1:NN-RN for j=1:1:NN-RN

ANa(i,j)=Ka(M(i),M(j)); bna(i,1)=pxa(M(i),1); end end

for i=1:RN

for j=1:NN-RN

BNa(i,j)=Ka(RC(i),M(j)); end end

AN=cell2mat(ANa); bn=cell2mat(bna); BN=cell2mat(BNa);

X=Jordan(AN,bn); %利用Gauss-Jordan法求解非约束结点的位移X %%----------------求解非约束结点的位移X-------------% %----------------由X可得所有结点位移x-------------% BN=BN*X; m=1; n=1; for i=1:1:NN if x(2*i)==1

x(2*i-1)=X(m); x(2*i)=X(m+1); m=m+2; else if x(2*i)==0

px(2*i-1)=BN(n); px(2*i)=BN(n+1); n=n+2; end

有限元课作业,我还是MATLAB稍微会点,C好久不用就没用C写,这个是大作业报告。

end end

%% 列主元Jordan消去法 将系数矩阵化成对角矩阵求解方程组的数值解 function x=Jordan(A,b) %开始计算，赋初值 [n,m]=size(A); x=zeros(n,1); for k=1:n

%选主元 max1=0; for i=k:n

if abs(A(i,k))>max1 max1=abs(A(i,k)); r=i; …… 此处隐藏：2951字，全部文档内容请下载后查看。喜欢就下载吧 ……