%用四節(jié)點平面板單元求解平板中心受集中載荷問題
%------------------------------------
%  輸入控制參數(shù)
%------------------------------------

clear                    %清除workplace殘留數(shù)據(jù)
nel=4;                   % 單元數(shù)
nnel=4;                  % 每個單元的節(jié)點數(shù)
ndof=3;                  % 每個節(jié)點的自由度數(shù)
nnode=9;                 % 系統(tǒng)總節(jié)點數(shù)
sdof=nnode*ndof;         % 系統(tǒng)總自由度數(shù)  
edof=nnel*ndof;          % 每個單元的自由度數(shù)
emodule=30e6;            % 楊氏彈性模量
poisson=0.3;             % 泊松比
t=0.1;                   % 薄板厚度
nglxb=2; nglyb=2;        % 彎曲對應(yīng)的2x2高斯拉格朗日積分 
nglb=nglxb*nglyb;        % 彎曲對應(yīng)的每個單元的高斯積分點數(shù)
nglxs=1; nglys=1;        % 剪切對應(yīng)的1x1高斯拉格朗日積分 
ngls=nglxs*nglys;        % 剪切對應(yīng)的每個單元的高斯積分點數(shù)

%---------------------------------------------
%  輸入節(jié)點坐標(biāo)值
%  gcoord(i,j) i：節(jié)點號　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];

%---------------------------------------------------------
%  每個單元對應(yīng)的節(jié)點號（逆時針排列）
%  nodes(i,j) i:節(jié)點號 j:對應(yīng)的單元號
%---------------------------------------------------------

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);        % 對應(yīng)的值

%----------------------------------------------
% 初始化矩陣和矢量
%----------------------------------------------

ff=zeros(sdof,1);       % 載荷矢量
kk=zeros(sdof,sdof);    % 系統(tǒng)剛度矩陣
disp=zeros(sdof,1);     % 系統(tǒng)位移矢量
index=zeros(edof,1);    % 每個單元所包含的自由度
kinmtpb=zeros(3,edof);   % 彎曲幾何函數(shù)矩陣
matmtpb=zeros(3,3);      % 彎曲材料系數(shù)矩陣
kinmtps=zeros(2,edof);   % 剪切幾何函數(shù)矩陣
matmtps=zeros(2,2);      % 剪切材料系數(shù)矩陣

%----------------------------
% 載荷矢量
%----------------------------

ff(27)=10;              % 結(jié)點9所受的集中載荷

%-----------------------------------------------------------------
%  單元剛度矩陣計算及其組合
%-----------------------------------------------------------------
%
%  彎曲相關(guān)計算
%
[pointb,weightb]=swp2(nglxb,nglyb);     % 積分點和權(quán)系數(shù)
matmtpb=sbm(emodule,poisson)*t^3/12;  %彎曲材料系數(shù)
%
%  剪切相關(guān)計算
%
[points,weights]=swp2(nglxs,nglys);     % 積分點和權(quán)系數(shù)
shearm=0.5*emodule/(1.0+poisson);          % 剪切模量
shcof=5/6;                                 % 剪切修正因數(shù) 
matmtps=shearm*shcof*t*[1 0; 0 1];         % 剪切材料系數(shù)矩陣

for iel=1:nel           % 對所有單元數(shù)的循環(huán)

for i=1:nnel
nd(i)=nodes(iel,i);         % 當(dāng)前單元對應(yīng)的節(jié)點
xcoord(i)=gcoord(nd(i),1);  % 節(jié)點對應(yīng)的x坐標(biāo)值
ycoord(i)=gcoord(nd(i),2);  % 節(jié)點對應(yīng)的y坐標(biāo)值
end

k=zeros(edof,edof);         % 初始化單元剛度矩陣
kb=zeros(edof,edof);        % 初始化彎曲剛度矩陣
ks=zeros(edof,edof);        % 初始化剪切剛度矩陣

%------------------------------------------------------
%  彎曲相關(guān)計算
%------------------------------------------------------

for intx=1:nglxb
x=pointb(intx,1);                  % x軸積分點坐標(biāo)
wtx=weightb(intx,1);               % 權(quán)系數(shù)
for inty=1:nglyb
y=pointb(inty,2);                  % y軸積分點坐標(biāo)
wty=weightb(inty,2) ;              % 權(quán)系數(shù)

[shape,dhdr,dhds]=ssf(x,y);     % 計算形函數(shù)和對其相應(yīng)的求導(dǎo)

jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord);  % 計算雅可比行列式

detjacob=det(jacob2);                 % 計算雅可比行列式的值
invjacob=inv(jacob2);                 % 求雅可比行列式的逆

[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 計算ddhdr,dhds在迪卡爾坐標(biāo)下的值

kinmtpb=sbB(nnel,dhdx,dhdy);          % 計算彎曲幾何函數(shù)矩陣

%--------------------------------------------
%  計算彎曲剛度矩陣
%--------------------------------------------

kb=kb+kinmtpb'*matmtpb*kinmtpb*wtx*wty*detjacob;

end
end                      % 結(jié)束彎曲剛度矩陣的計算

%------------------------------------------------------
% 剪切相關(guān)計算
%------------------------------------------------------

for intx=1:nglxs
x=points(intx,1);                  % x軸積分點坐標(biāo)
wtx=weights(intx,1);               % 權(quán)系數(shù)
for inty=1:nglys
y=points(inty,2);                  % y軸積分點坐標(biāo)
wty=weights(inty,2) ;              % 權(quán)系數(shù)

[shape,dhdr,dhds]=ssf(x,y);     % 計算形函數(shù)和對其相應(yīng)的教學(xué)求導(dǎo)

jacob2=sjc(nnel,dhdr,dhds,xcoord,ycoord);  % 計算雅可比行列式

detjacob=det(jacob2);                 % 計算雅可比行列式的值
invjacob=inv(jacob2);                 % 求雅可比行列式的逆

[dhdx,dhdy]=sxy(nnel,dhdr,dhds,invjacob); % 計算dhdr,dhds在迪卡爾坐標(biāo)下的值

kinmtps=ssB(nnel,dhdx,dhdy,shape);        % 計算剪切幾何函數(shù)矩陣

%----------------------------------------
%  計算剪切剛度矩陣
%----------------------------------------

ks=ks+kinmtps'*matmtps*kinmtps*wtx*wty*detjacob;   

end
end                      % 結(jié)束剪切剛度矩陣的計算

%--------------------------------
% 計算單元  剛度矩陣
%--------------------------------

k=kb+ks;

index=etsd(nd,nnel,ndof);% 單元對應(yīng)的系統(tǒng)自由度號

kk=ask(kk,k,index);  % 合成系統(tǒng)剛度矩陣

end

%-----------------------------
%   加載邊界條件
%-----------------------------

[kk,ff]=dbc(kk,ff,bcdof,bcval);

%----------------------------
% 求解
%----------------------------

disp=kk\ff;   

num=1:1:sdof;
nodedisp=[num' disp]                  % 輸出節(jié)點位移
%----------------------------
% 后處理
%----------------------------
result=zeros(25,3);                   %初始化
displac(75)=0;　　　　　　　　　　　 　%將輸出的結(jié)果從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);             %將節(jié)點位移以節(jié)點順序輸出，i->節(jié)點號，j->對應(yīng)位移
exgcoord=excoord;　　　　　　　　　　　%將節(jié)點坐標(biāo)從5x5的四分之一板擴充到10x10的全板
[finresult]=agr(exgcoord,result);　　 %將全板的節(jié)點坐標(biāo)和節(jié)點位移對應(yīng)起來，i->節(jié)點號，j->對應(yīng)坐標(biāo)和位移
for i=1:25;
    finresultin(1,i)=sqrt(finresult(3,i)^2+finresult(4,i)^2+finresult(5,i)^2);　　%求節(jié)點位移
end
Z=arrayfin(finresultin);　　　　　　　%排列節(jié)點位移
surf(Z);　　　　　　　　　　　　　　　　%畫板變形圖