8节点等参单元matlab编程
function u=u() % 本程序采用8结点四边形单元，分析固定位移边界条件的弹性力学平面问题fprintf('请输入弹性模量E的值(单位均采用国际单位制，下同) \n ');
E=input('');
fprintf('请输入泊松比NU的值\n ');
NU=input('');
fprintf('请输入结构的厚度hou\n ');
hou=input('');
fprintf('请输入分析类型p，平面应力问题请输入1，平面应变问题请输入2\n ');
p=input('');
while p~=1&p~=2
fprintf('输入有误。请重新输入分析类型p，平面应力问题请输入1，平面应变问题请输入2\n ');
p=1;
end
if p==1 % 根据问题类型确定弹性矩阵D
D=(E/(1-NU*NU))*;
elseif p==2
D=(E/(1+NU)/(1-2*NU))*;
end
fprintf('请输入结点个数nn(nn为不小于8的正整数) \n ');
nn=input('');
while nn<8
fprintf('输入有误，请重新输入结点个数nn(nn为不小于8的正整数) \n ');
nn=input('');
end
fprintf('请输入单元个数ne(ne为不小于1的正整数) \n ');
ne=input('');
while ne<1
fprintf('输入有误，请重新输入结点个数ne(ne为不小于1的正整数) \n ');
ne=input('');
end
K1=zeros(2*nn); % K1为结构整体刚度矩阵
for i=1:ne
fprintf('请输入第%d个单元的结点坐标co(形式为);\n',i);
co=input('');
fprintf('请输入该单元的结点编号no(形式为) \n ');
no=input('');
k1=k(co,E,NU,D,hou); % 调用单元刚度矩阵函数
K1=kk(no,k1,K1); % 将单元刚度矩阵装配到整体刚度矩阵中
end
fprintf('未修正位移边界条件的整体刚度矩阵为K1=\n');
disp(K1);
fprintf('请输入各结点外荷载值(形式为,约束端的外荷载按0输入) \n ');
F=input('');
while size(F)~=2*nn
fprintf('输入有误。请重新输入各结点外荷载值(形式为) \n ');
F=input('');
end
%以下采用对角元素乘大数法修正位移边界条件
fprintf('请输入x方向受到位移约束的结点编号及其约束位移值(形式为) \n ');
xd=input('');
xn=size(xd,1);
AA=1.0e30; % AA为大数
for i=1:xn
K1(2*xd(i,1)-1,2*xd(i,1)-1)=K1(2*xd(i,1)-1,2*xd(i,1)-1)*AA;
F(2*xd(i,1)-1)=xd(i,2)*K1(2*xd(i,1)-1,2*xd(i,1)-1);
end
fprintf('请输入y方向受到位移约束的结点编号及其约束位移值(形式为) \n ');
yd=input('');
yn=size(yd,1);
for i=1:yn
K1(2*yd(i,1),2*yd(i,1))=K1(2*yd(i,1),2*yd(i,1))*AA;
F(2*yd(i,1))=yd(i,2)*K1(2*yd(i,1),2*yd(i,1));
end
u=K1\F; % 解整体刚度方程KU=P
fprintf('边界修正后的整体刚度矩阵为K1=\n');
disp(K1);
fprintf('边界修正后的外荷载列为F=\n');
disp(F);
fprintf('在外荷载作用下各结点的位移为u=\n');
disp(u);
〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓
function k=k(co,E,NU,D,hou)
% co为单元的结点坐标，为8行2列矩阵
% 函数返回单元的刚度矩阵
syms s t; % s,t表示局部坐标
N=sym(zeros(8,1));
N(1)=(1-s)*(1-t)*(-s-t-1)/4;
N(2)=(1+s)*(1-t)*(s-t-1)/4;
N(3)=(1+s)*(1+t)*(s+t-1)/4;
N(4)=(1-s)*(1+t)*(-s+t-1)/4;
N(5)=(1-t)*(1+s)*(1-s)/2;
N(6)=(1+s)*(1+t)*(1-t)/2;
N(7)=(1+t)*(1+s)*(1-s)/2;
N(8)=(1-s)*(1+t)*(1-t)/2;
x=0;
y=0;
for i=1:8
x=x+N(i)*co(i,1);
y=y+N(i)*co(i,2);
end
xs=diff(x,s);
xt=diff(x,t);
ys=diff(y,s);
yt=diff(y,t);
J=xs*yt-ys*xt;
for i=1:8
Ns(i)=diff(N(i),s);
Nt(i)=diff(N(i),t);
end
for i=1:8
g(i)=yt*Ns(i)-ys*Nt(i);
h(i)=xs*Nt(i)-xt*Ns(i);
end
B=sym(zeros(3,16));
for i=1:8
B(1,2*i-1)=B(1,2*i-1)+g(i);
B(2,2*i)=B(2,2*i)+h(i);
B(3,2*i-1)=B(3,2*i-1)+h(i);
B(3,2*i)=B(3,2*i)+g(i);
end
Bnew=simplify(B);
Jnew=simplify(J);
BD=transpose(Bnew)*D*Bnew/Jnew;
r=int(int(BD,t,-1,1),s,-1,1);
z=hou*r;
k=double(z);
〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓
function K=kk(no,k,K)
%no表示相应单元在整体坐标系中的编号，k为该单元的刚度矩阵
%K返回结构的整体刚度矩阵
for i=1:8
for j=1:8
K(2*no(i)-1,2*no(j)-1)=K(2*no(i)-1,2*no(j)-1)+k(2*i-1,2*j-1);
K(2*no(i)-1,2*no(j))=K(2*no(i)-1,2*no(j))+k(2*i-1,2*j);
K(2*no(i),2*no(j)-1)=K(2*no(i),2*no(j)-1)+k(2*i,2*j-1);
K(2*no(i),2*no(j))=K(2*no(i),2*no(j))+k(2*i,2*j);
end
end
〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓
function u=u() % 本程序采用8结点四边形单元，分析固定位移边界条件的弹性力学平面问题
fprintf('请输入弹性模量E的值(单位均采用国际单位制，下同) \n ');
E=input('');
fprintf('请输入泊松比NU的值\n ');
NU=input('');
fprintf('请输入结构的厚度hou\n ');
hou=input('');
fprintf('请输入分析类型p，平面应力问题请输入1，平面应变问题请输入2\n ');
p=input('');
while p~=1&p~=2
fprintf('输入有误。请重新输入分析类型p，平面应力问题请输入1，平面应变问题请输入2\n ');
p=1;
end
if p==1 % 根据问题类型确定弹性矩阵D
D=(E/(1-NU*NU))*;
elseif p==2
D=(E/(1+NU)/(1-2*NU))*;
end
fprintf('请输入结点个数nn(nn为不小于8的正整数) \n ');
nn=input('');
while nn<8
fprintf('输入有误，请重新输入结点个数nn(nn为不小于8的正整数) \n ');
nn=input('');
end
fprintf('请输入单元个数ne(ne为不小于1的正整数) \n ');
ne=input('');
while ne<1
fprintf('输入有误，请重新输入结点个数ne(ne为不小于1的正整数) \n ');
ne=input('');
end
K1=zeros(2*nn); % K1为结构整体刚度矩阵
for i=1:ne
fprintf('请输入第%d个单元的结点坐标co(形式为);\n',i);
co=input('');
fprintf('请输入该单元的结点编号no(形式为) \n ');
no=input('');
k1=k(co,E,NU,D,hou); % 调用单元刚度矩阵函数
K1=kk(no,k1,K1); % 将单元刚度矩阵装配到整体刚度矩阵中
end
fprintf('未修正位移边界条件的整体刚度矩阵为K1=\n');
disp(K1);
fprintf('请输入各结点外荷载值(形式为,约束端的外荷载按0输入) \n ');
F=input('');
while size(F)~=2*nn
fprintf('输入有误。请重新输入各结点外荷载值(形式为) \n ');
F=input('');
end
%以下采用对角元素乘大数法修正位移边界条件
fprintf('请输入x方向受到位移约束的结点编号及其约束位移值(形式为) \n ');
xd=input('');
xn=size(xd,1);
AA=1.0e30; % AA为大数
for i=1:xn
K1(2*xd(i,1)-1,2*xd(i,1)-1)=K1(2*xd(i,1)-1,2*xd(i,1)-1)*AA;
F(2*xd(i,1)-1)=xd(i,2)*K1(2*xd(i,1)-1,2*xd(i,1)-1);
end
fprintf('请输入y方向受到位移约束的结点编号及其约束位移值(形式为) \n ');
yd=input('');
yn=size(yd,1);
for i=1:yn
K1(2*yd(i,1),2*yd(i,1))=K1(2*yd(i,1),2*yd(i,1))*AA;
F(2*yd(i,1))=yd(i,2)*K1(2*yd(i,1),2*yd(i,1));
end
u=K1\F; % 解整体刚度方程KU=P
fprintf('边界修正后的整体刚度矩阵为K1=\n');
disp(K1);
fprintf('边界修正后的外荷载列为F=\n');
disp(F);
fprintf('在外荷载作用下各结点的位移为u=\n');
disp(u);
〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓
function k=k(co,E,NU,D,hou)
% co为单元的结点坐标，为8行2列矩阵
% 函数返回单元的刚度矩阵
syms s t; % s,t表示局部坐标
N=sym(zeros(8,1));
N(1)=(1-s)*(1-t)*(-s-t-1)/4;
N(2)=(1+s)*(1-t)*(s-t-1)/4;
N(3)=(1+s)*(1+t)*(s+t-1)/4;
N(4)=(1-s)*(1+t)*(-s+t-1)/4;
N(5)=(1-t)*(1+s)*(1-s)/2;
N(6)=(1+s)*(1+t)*(1-t)/2;
N(7)=(1+t)*(1+s)*(1-s)/2;
N(8)=(1-s)*(1+t)*(1-t)/2;
x=0;
y=0;
for i=1:8
x=x+N(i)*co(i,1);
y=y+N(i)*co(i,2);
end
xs=diff(x,s);
xt=diff(x,t);
ys=diff(y,s);
yt=diff(y,t);
J=xs*yt-ys*xt;
for i=1:8
Ns(i)=diff(N(i),s);
Nt(i)=diff(N(i),t);
end
for i=1:8
g(i)=yt*Ns(i)-ys*Nt(i);
h(i)=xs*Nt(i)-xt*Ns(i);
end
B=sym(zeros(3,16));
for i=1:8
B(1,2*i-1)=B(1,2*i-1)+g(i);
B(2,2*i)=B(2,2*i)+h(i);
B(3,2*i-1)=B(3,2*i-1)+h(i);
B(3,2*i)=B(3,2*i)+g(i);
end
Bnew=simplify(B);
Jnew=simplify(J);
BD=transpose(Bnew)*D*Bnew/Jnew;
r=int(int(BD,t,-1,1),s,-1,1);
z=hou*r;
k=double(z);
〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓〓
function K=kk(no,k,K)
%no表示相应单元在整体坐标系中的编号，k为该单元的刚度矩阵
%K返回结构的整体刚度矩阵
for i=1:8
for j=1:8
K(2*no(i)-1,2*no(j)-1)=K(2*no(i)-1,2*no(j)-1)+k(2*i-1,2*j-1);
K(2*no(i)-1,2*no(j))=K(2*no(i)-1,2*no(j))+k(2*i-1,2*j);
K(2*no(i),2*no(j)-1)=K(2*no(i),2*no(j)-1)+k(2*i,2*j-1);
K(2*no(i),2*no(j))=K(2*no(i),2*no(j))+k(2*i,2*j);
end
end 写得不错 先下下来 需要用到 赞一个 看来比较简单,不知道写的怎么样啊? 恩，谢谢楼主的无私分享，学习一下
页:
[1]