设为首页收藏本站联系我们
切换到宽版

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛»论坛 › 算法编程 › Matlab讨论区 › 8节点等参单元matlab编程
返回列表 发新帖
查看: 2090|回复: 4

[共享资源] 8节点等参单元matlab编程

[复制链接]
huweifa
发表于 2006-8-12 18:47 | 显示全部楼层 |阅读模式

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要 登录 才可以下载或查看，没有账号？我要加入

x
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))*[1 NU 0; NU 1 0;0 0 (1-NU)/2];
elseif p==2
   D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
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(形式为[x1 y1;x2 y2;……;x8 y8]);\n',i);
    co=input('');
    fprintf('请输入该单元的结点编号no(形式为[1 2……7 8]) \n ');
    no=input('');
    k1=k(co,E,NU,D,hou);          % 调用单元刚度矩阵函数
    K1=kk(no,k1,K1);              % 将单元刚度矩阵装配到整体刚度矩阵中
end

fprintf('未修正位移边界条件的整体刚度矩阵为K1=\n');
disp(K1);

fprintf('请输入各结点外荷载值(形式为[fx1;fy1;fx2;fy2;……;fxn;fyn],约束端的外荷载按0输入) \n ');
F=input('');
while size(F)~=2*nn
    fprintf('输入有误。请重新输入各结点外荷载值(形式为[fx1;fy1;fx2;fy2;……;fxn;fyn]) \n ');
    F=input('');
end
   
                                  %  以下采用对角元素乘大数法修正位移边界条件
fprintf('请输入x方向受到位移约束的结点编号及其约束位移值(形式为[node1 u1;node2 u2;……;noden un]) \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方向受到位移约束的结点编号及其约束位移值(形式为[node1 v1;node2 v2;……;noden vn]) \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))*[1 NU 0; NU 1 0;0 0 (1-NU)/2];
elseif p==2
   D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
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(形式为[x1 y1;x2 y2;……;x8 y8]);\n',i);
    co=input('');
    fprintf('请输入该单元的结点编号no(形式为[1 2……7 8]) \n ');
    no=input('');
    k1=k(co,E,NU,D,hou);          % 调用单元刚度矩阵函数
    K1=kk(no,k1,K1);              % 将单元刚度矩阵装配到整体刚度矩阵中
end

fprintf('未修正位移边界条件的整体刚度矩阵为K1=\n');
disp(K1);

fprintf('请输入各结点外荷载值(形式为[fx1;fy1;fx2;fy2;……;fxn;fyn],约束端的外荷载按0输入) \n ');
F=input('');
while size(F)~=2*nn
    fprintf('输入有误。请重新输入各结点外荷载值(形式为[fx1;fy1;fx2;fy2;……;fxn;fyn]) \n ');
    F=input('');
end
   
                                  %  以下采用对角元素乘大数法修正位移边界条件
fprintf('请输入x方向受到位移约束的结点编号及其约束位移值(形式为[node1 u1;node2 u2;……;noden un]) \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方向受到位移约束的结点编号及其约束位移值(形式为[node1 v1;node2 v2;……;noden vn]) \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

查看全部评分

回复
分享到：

使用道具 举报

argye
发表于 2006-8-12 19:57 | 显示全部楼层
写得不错
回复 支持 反对

使用道具 举报

Gaviny
发表于 2006-9-12 22:13 | 显示全部楼层
先下下来   需要用到    赞一个
回复 支持 反对

使用道具 举报

fx2003521
发表于 2006-9-14 11:51 | 显示全部楼层
看来比较简单,不知道写的怎么样啊?
回复 支持 反对

使用道具 举报

星矢
发表于 2011-7-20 13:04 | 显示全部楼层
恩，谢谢楼主的无私分享，学习一下
回复 支持 反对

使用道具 举报

返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-11 05:19 , Processed in 0.098553 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表