转：我的有限元学习笔记 (Part 2 - 用Matlab编写的DSM方法范例程序)

我的有限元学习笔记 (Part 2 - 用Matlab编写的DSM方法范例程序)

[ 2004-07-18 18:20:26 | Author: moneywood ]

Font Size: Large | Medium | Small

本文是我的有限元学习笔记的一部分，将随我学习进度在blog陆续po出，本人新手，请各位看官多提意见，谢谢！

Part 2 ...
这段程序是根据第五章中提供的Mathematica程序改写的，说实话Mathematica对于500行内的小规模试算和验证效率大大超过Matlab，不过我还是做了这个么个改写，目的就是练手，并比较两个软件在运作方式上的异同。
当然，我改写了以后，这个程序只能进行数值计算了，这个没办法，不过我在本来程序的基础上做了一些功能的添加，我觉得更加 general,user-friendly了，在主程序部分，你只要按顺序输入连接点的坐标(2D)，以及连接后Element属性(连接点，截面 积)，主刚度阵的生成非常简便，不过BC你需要自己定义，毕竟这么个小程序，只能做很简单的事儿，要求别太高了:)

function ifem( input_args )
%IFEM Summary of this function goes here
% DSM of Example Truss
% Trnsferd from Mathematica
Km=zeros(6); % Primary Stiffness Matrix
P=[0,0;10,0;10,10]; % point coordinate
% Element proporties of each element
E1=[1,3,2*sqrt(2)];
E2=[1,2,1];
E3=[2,3,0.5];
% generate primary stiffness matrix
Ke1=ElemStiff2DTwoNodeBar(E1,P);
Km=MergeElemIntoMasterStiff(Ke1, E1, Km);
Ke2=ElemStiff2DTwoNodeBar(E2,P);
Km=MergeElemIntoMasterStiff(Ke2, E2, Km);
Ke3=ElemStiff2DTwoNodeBar(E3,P);
Km=MergeElemIntoMasterStiff(Ke3, E3, Km);
disp('Master Stiff Matrix :');disp(Km);

% apply BC
Kmod=ModifiedMasterStiffForDBC([1,2,4],Km);
disp('Modified Master Stiff Matrix :');disp(Kmod);
f=([0 0 0 0 2 1])';
% solute
u=inv(Kmod)*f;
disp('displacement : ');disp(u');

% recover internal force
disp(sprintf(' Internal Force of Element 1 : %f', IntForce2DTwoNodeBar(E1, u, P)));
disp(sprintf(' Internal Force of Element 2 : %f', IntForce2DTwoNodeBar(E2, u, P)));
disp(sprintf(' Internal Force of Element 3 : %f', IntForce2DTwoNodeBar(E3, u, P)));


%-----------------------------------------------------
function Ke = ElemStiff2DTwoNodeBar( Elem, P )
Em=100;
dx=P(Elem(2),1)-P(Elem(1),1);
dy=P(Elem(2),2)-P(Elem(1),2);
L=sqrt(dx^2+dy^2);
c=dx/L;
s=dy/L;
Ke=[c^2, c*s,-c^2,-c*s;
c*s, s^2,-s*c,-s^2;
-c^2,-s*c, c^2, s*c;
-s*c,-s^2, s*c, s^2];
Ke=Ke*(Em*Elem(3)/L);
end

function Km = MergeElemIntoMasterStiff( Ke, Elem, Km )
eftab=[Elem(1)*2-1,Elem(1)*2,Elem(2)*2-1,Elem(2)*2];
for i=1:4
ii=eftab(i);
for j=i:4
jj=eftab(j);
Km(jj,ii) = Km(jj,ii) + Ke(i,j);
Km(ii,jj) = Km(jj,ii);
end
end

function Kmod=ModifiedMasterStiffForDBC (pdof, Km)
Kmod=Km;
for k=1:1:length(pdof)
i=pdof(k);
Kmod(i,1:length(Kmod))=0;
Kmod(1:length(Kmod),i)=0;
Kmod(i,i)=1;
end

function fmod=ModifiedMasterForcesForDBC (pdof, f)
fmod=f;
for k=1:1:length(pdof)
i=pdof(k);
fmod(i)=0;
end

function p=IntForce2DTwoNodeBar(Elem, u, P)
Em=100;
dx=P(Elem(2),1)-P(Elem(1),1);
dy=P(Elem(2),2)-P(Elem(1),2);
L=sqrt(dx^2+dy^2);
c=dx/L;
s=dy/L;
edtab=[Elem(1)*2-1,Elem(1)*2,Elem(2)*2-1,Elem(2)*2];
ubar=[c*u(edtab(1))+s*u(edtab(2)),-s*u(edtab(1))+c*u(edtab(2)),c*u(edtab(3))+s*u(edtab(4)),-s*u(edtab(3))+c*u(edtab(4))];
e=(ubar(3)-ubar(1))/L;
p=Em*Elem(3)*e;