8节点等参单元matlab编程.txt

本程序采用8结点四边形单元

8节点等参单元matlab编程 
 
作者：huweifa    文章来源：本站原创    点击数：425    更新时间：2006-8-12 
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