%用四节点平面板单元求解平板中心受集中载荷问题
%------------------------------------
% 输入控制参数
%------------------------------------
clear %清除workplace残留数据
nel=4; % 单元数
nnel=4; % 每个单元的节点数
ndof=3; % 每个节点的自由度数
nnode=9; % 系统总节点数
sdof=nnode*ndof; % 系统总自由度数
edof=nnel*ndof; % 每个单元的自由度数
emodule=30e6; % 杨氏弹性模量
poisson=0.3; % 泊松比
t=0.1; % 薄板厚度
nglxb=2; nglyb=2; % 弯曲对应的2x2高斯拉格朗日积分
nglb=nglxb*nglyb; % 弯曲对应的每个单元的高斯积分点数
nglxs=1; nglys=1; % 剪切对应的1x1高斯拉格朗日积分
ngls=nglxs*nglys; % 剪切对应的每个单元的高斯积分点数
%---------------------------------------------
% 输入节点坐标值
% gcoord(i,j) i：节点号 j：x,y值
%---------------------------------------------
gcoord=[0.0 0.0; 2.5 0.0; 5.0 0.0;
0.0 2.5; 2.5 2.5; 5.0 2.5;
0.0 5.0; 2.5 5.0; 5.0 5.0];
%---------------------------------------------------------
% 每个单元对应的节点号（逆时针排列）
% nodes(i,j) i:节点号 j:对应的单元号
%---------------------------------------------------------
nodes=[1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8];
%-------------------------------------
% 输入边界条件
%-------------------------------------
bcdof=[1 2 3 4 6 7 9 11 12 16 20 21 23 25 26]; % 约束的自由度
bcval=zeros(1,15); % 对应的值
%----------------------------------------------
% 初始化矩阵和矢量
%----------------------------------------------
ff=zeros(sdof,1); % 载荷矢量
kk=zeros(sdof,sdof); % 系统刚度矩阵
disp=zeros(sdof,1); % 系统位移矢量
index=zeros(edof,1); % 每个单元所包含的自由度
kinmtpb=zeros(3,edof); % 弯曲几何函数矩阵
matmtpb=zeros(3,3); % 弯曲材料系数矩阵
kinmtps=zeros(2,edof); % 剪切几何函数矩阵
matmtps=zeros(2,2); % 剪切材料系数矩阵
%----------------------------
% 载荷矢量
%----------------------------
ff(27)=10; % 结点9所受的集中载荷
%-----------------------------------------------------------------
% 单元刚度矩阵计算及其组合
%-----------------------------------------------------------------
%
% 弯曲相关计算
%
[pointb,weightb]=swp2(nglxb,nglyb); % 积分点和权系数
matmtpb=sbm(emodule,poisson)*t^3/12; %弯曲材料系数
%
% 剪切相关计算
%
[points,weights]=swp2(nglxs,nglys); % 积分点和权系数
shearm=0.5*emodule/(1.0+poisson); % 剪切模量
shcof=5/6; % 剪切修正因数
matmtps=shearm*shcof*t*[1 0; 0 1]; % 剪切材料系数矩阵
for iel=1:nel % 对所有单元数的循环
for i=1:nnel
nd(i)=nodes(iel,i); % 当前单元对应的节点
xcoord(i)=gcoord(nd(i),1); % 节点对应的x坐标值
ycoord(i)=gcoord(nd(i),2); % 节点对应的y坐标值
end
k=zeros(edof,edof); % 初始化单元刚度矩阵
kb=zeros(edof,edof); % 初始化弯曲刚度矩阵
ks=zeros(edof,edof); % 初始化剪切刚度矩阵
%------------------------------------------------------
% 弯曲相关计算
%------------------------------------------------------
for intx=1:nglxb
x=pointb(intx,1); % x轴积分点坐标
wtx=weightb(intx,1); % 权系数
for inty=1:nglyb
y=pointb(inty,2); % y轴积分点坐标
wty=weightb(inty,2) ; % 权系数
[shape,dhdr,dhds]=ssf(x,y); % 计算形函数和对其相应的求导
jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord); % 计算雅可比行列式
detjacob=det(jacob2); % 计算雅可比行列式的值
invjacob=inv(jacob2); % 求雅可比行列式的逆
[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 计算ddhdr,dhds在迪卡尔坐标下的值
kinmtpb=sbB(nnel,dhdx,dhdy); % 计算弯曲几何函数矩阵
%--------------------------------------------
% 计算弯曲刚度矩阵
%--------------------------------------------
kb=kb+kinmtpb'*matmtpb*kinmtpb*wtx*wty*detjacob;
end
end % 结束弯曲刚度矩阵的计算
%------------------------------------------------------
% 剪切相关计算
%------------------------------------------------------
for intx=1:nglxs
x=points(intx,1); % x轴积分点坐标
wtx=weights(intx,1); % 权系数
for inty=1:nglys
y=points(inty,2); % y轴积分点坐标
wty=weights(inty,2) ; % 权系数
[shape,dhdr,dhds]=ssf(x,y); % 计算形函数和对其相应的教学求导
jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord); % 计算雅可比行列式
detjacob=det(jacob2); % 计算雅可比行列式的值
invjacob=inv(jacob2); % 求雅可比行列式的逆
[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 计算dhdr,dhds在迪卡尔坐标下的值
kinmtps=ssB(nnel,dhdx,dhdy,shape); % 计算剪切几何函数矩阵
%----------------------------------------
% 计算剪切刚度矩阵
%----------------------------------------
ks=ks+kinmtps'*matmtps*kinmtps*wtx*wty*detjacob;
end
end % 结束剪切刚度矩阵的计算
%--------------------------------
% 计算单元 刚度矩阵
%--------------------------------
k=kb+ks;
index=etsd(nd,nnel,ndof);% 单元对应的系统自由度号
kk=ask(kk,k,index); % 合成系统刚度矩阵
end
%-----------------------------
% 加载边界条件
%-----------------------------
[kk,ff]=dbc(kk,ff,bcdof,bcval);
%----------------------------
% 求解
%----------------------------
disp=kk\ff;
num=1:1:sdof;
nodedisp=[num' disp] % 输出节点位移
%----------------------------
% 后处理
%----------------------------
result=zeros(25,3);
displace(75)=0;
a=re(1,0,disp);
a(75)=0;
displace=displace+a;
a=re(10,6,disp);
a(75)=0;
displace=displace+a;
a=re(19,12,disp);
a(75)=0;
displace=displace+a;
for i=1:15;
displace(45+i)=displace(15+i);
displace(60+i)=displace(i);
end
[result]=dtxy(displace);
for i=1:25
displacex(i,:)=result(:,1);
displacex(i,:)=result(:,1);
displacey(i,:)=result(:,2);
displacez(i,:)=result(:,3);
end
[X,Y] = meshgrid(0:10/24:10);
Z=sqrt(displacex.^2+displacey.^2+displacez.^2);
surf(Z);
