平面四节点等参单元matlab实现(4)

约束的节点

n=size(X,1); %节点数 m=size(elem,1); %单元数

K=zeros(2*n); %初始总体刚度矩阵 for i=1:m

syms Ks Et x y I1 I2 a b B; %定义可能存在的变量 e=[1,1;-1,1;-1,-1;1,-1];

for j=1:4 %形函数 N(j)=0.25*(1+e(j,1)*Ks)*(1+e(j,2)*Et); end x=0;y=0;

for j=1:4 %标准母单元映射到真实单元 x=x+N(j)*X(elem(i,j),1); y=y+N(j)*X(elem(i,j),2); end

J1=jacobian([x;y],[Ks;Et]); %雅克比矩阵及其转置 J=J1';

for j=1:4

I1=diff(N(j),Ks); %形函数分别对Ks和Et的偏导数 I2=diff(N(j),Et); C=(J^-1)*[I1;I2];

a=C(1); %形函数对x，y的偏导数 b=C(2);

B(1,2*j-1)=a; %组成B阵 B(1,2*j)=0; B(2,2*j-1)=0; B(2,2*j)=b; B(3,2*j-1)=b; B(3,2*j)=a; end

D=(E/(1-NU*NU))*[1,NU,0;NU,1,0;0,0,(1-NU)/2]; %D阵 k=zeros(8,8);

Kss=[-0.906179,-0.538469,0,0.538469,0.906179]; %5*5高斯积分点 Ett=[-0.906179,-0.538469,0,0.538469,0.906179];

H=[0.236926,0.478628,0.568888,0.478628,0.236926];%高斯积分权系数 for j=1:5 %高斯积分求单元刚

度阵

for l=1:5

Ks=Kss(j);Et=Ett(l); B=subs(B); J=subs(J);