untitled5.m

力学计算 力学计算 力学计算 力学计算 力学计算

%用四节点平面板单元求解平板中心受集中载荷问题
%------------------------------------
%  输入控制参数
%------------------------------------

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);                   %初始化
displac(75)=0;　　　　　　　　　　　 　%将输出的结果从5x5的四分之一板扩充到10x10的全板
a=re(1,0,disp);
a(75)=0;
displac=displace+a;
a=re(10,6,disp);
a(75)=0;
displac=displace+a;
a=re(19,12,disp);
a(75)=0;
displac=displace+a;
for i=1:15;
    displac(45+i)=displac(15+i);
    displac(60+i)=displac(i);
end
[result]=dtxy(displac);             %将节点位移以节点顺序输出，i->节点号，j->对应位移
exgcoord=excoord;　　　　　　　　　　　%将节点坐标从5x5的四分之一板扩充到10x10的全板
[finresult]=agr(exgcoord,result);　　 %将全板的节点坐标和节点位移对应起来，i->节点号，j->对应坐标和位移
for i=1:25;
    finresultin(1,i)=sqrt(finresult(3,i)^2+finresult(4,i)^2+finresult(5,i)^2);　　%求节点位移
end
Z=arrayfin(finresultin);　　　　　　　%排列节点位移
surf(Z);　　　　　　　　　　　　　　　　%画板变形图