plane4.m

四节点矩形单元matlab有限元程序

%function plane stress
% This function contains gives out results for  such a structure having  
%such constant as  the moduius of elasticity E, Possion's ratio NU, 
%thickness t,the numble of elements NE,the numble of freedom of a element
%NF, the numble of all nodes NP,the numble of nodes of each element ND,the
%numble of restriction NR ,the numble of freedom of a elementNDF, the
%numble of total nodes N.

%%开始输入材料参数%%
E=210e6;
NU=0.300;
t=0.025;
p=1;
%%输入完毕%%

%%输入板信息：NF——结点自由度数 　 ND——单元结点数　  LENGTH——板长度　 BREADTH——板宽度%%
ND=4; 
NF=2;
LENGTH=0.50;
BREADTH=0.25;
%%输入完毕%%

%%输入网格划分信息：NEI——横向数   NEJ——纵向数%%
NEI=8;
NEJ=8;
%%输入完毕%%

NE=NEI*NEJ;  %%单元总数
NDF=ND*NF;   %%每个单元的总自由度数
NP=(NEI+1)*(NEJ+1);  %%节点总数
N=NF*NP;     %%自由度总数
[X,Y]=BilinearQuadElementcoordinate(NEI,NEJ,LENGTH,BREADTH);
ME=BilinearQuadElementME(ND,NEI,NEJ);
hold on;
BilinearQuadElementShowmesh(X,Y,ME,NE,ND);

%%%%%%%%%%%%%%%%%%开始输入约束条件
for I=1:N
    U(I)=888;
    P(I)=888;
end
U=BilinearQuadElementRestictU(NF,U,X,Y);
P=BilinearQuadElementRestrictP(NF,P,X,Y); 
%%约束条件输入完毕%%
%%输入信息结束%%




K=zeros(N,N);
for I=1:NE
    IS=BilinearQuadElementIS(ND,NF,I,ME);
    k=BilinearQuadElementStiffness(E,NU,t,X(ME(1,I)),Y(ME(1,I)), X(ME(2,I)),Y(ME(2,I)), X(ME(3,I)),Y(ME(3,I)),X(ME(4,I)),Y(ME(4,I)),p);
   % K=BilinearQuadAssemble(K,k,IS,NDF);
    K=BilinearQuadAssemble(K,k,ME(1,I),ME(2,I),ME(3,I),ME(4,I));
end
[P,U]=BilinearQuadElementCaculate(P,U,K,N);
for I=1:NP
    X(I)=X(I)+U(2*I-1);
    Y(I)=Y(I)+U(2*I);
end
%hold on;
%BilinearQuadElementShowDistortion(X,Y,ME,NE,ND);
%%求单元应力，主应力%%
%%单元数NE，单元各个节点位移u%%
STRESS=[];
for I=1:NE
    V=[U(2*ME(1,I)-1) U(2*ME(1,I)) U(2*ME(2,I)-1) U(2*ME(2,I)) U(2*ME(3,I)-1) U(2*ME(3,I)) U(2*ME(4,I)-1) U(2*ME(4,I))];
    u=V';
    SIGMA=BilinearQuadElementStresses(E,NU,X(ME(1,I)),Y(ME(1,I)),X(ME(2,I)),Y(ME(2,I)),X(ME(3,I)),Y(ME(3,I)),X(ME(4,I)),Y(ME(4,I)),p,u);
    STRESS1=BilinearQuadElementPStresses(SIGMA);
    STRESS=[STRESS,STRESS1];
end