p080123d3.m

有限元平面问题的求解相关程序。 包括质量矩阵

%p080123.m  用三胞体插值。
%将位移、频率结果进行比较，完备的程序，验证正确！
%——————基本变量——————
clear;clc;
E1=200e3;E2=E1;     %第一种情况：均匀材料，各向同性。
%E1=200e3;E2=E1/50; %第二种情况：最角上硬，中间软
%E2=200e3;E1=E2/50; %第三种情况：最角上软，中间硬
t=1;NU=0.3;p=1;                             %厚度、泊松比、平面问题的种类p，1为平面应力，2为应变
density1=5.0e-6;density2=5.0e-6;fj=100;     %密度，外载荷，全部单位都与克、毫米相适应
line=9;ar=3; %原始有限元的单元行数，列数，
fixnod=line+1;          %固定节点数目
freenod=(line+1)*ar;    %加边界条件后总自由节点数,其数值和边界条件密切相关
fstart=line+2;          %第一个自由节点号
%---------单胞组成-------
dsubline=3;dsubar=3;    %单胞中包含小单元的行列数
vf=4/9;vb=5/9;          %两种材料的体积比
b=10/dsubar;  %b为胞体中小单元边长，由结构决定，
%-----------------等效所使用的变量----------------
Ee=vf*E1+vb*E2;Er=1/(vf/E1+vb/E2);  %刚度均匀化法、柔度均匀化法得到的等效刚度
subline=3;subar=1;a=dsubar*b;        %等效后，等效单元行、列数、边长
nod=dsubline*subline+1;         %等效前，加有载荷的节点数目
subnod=subline+1;               %等效后，加有载荷的等效节点数目
subfreenod=(subline+1)*subar;   %加边界条件并等效后总自由节点数,其数值和边界条件密切相关
subfstart=subline+2;    %加边界条件并等效后第一个自由节点号
%-----------------------------------------------------
%-----------------------------------------------------
k1=BilinearQuadElementStiffness(E1,NU,t,0,0,b,0,b,b,0,b,p);%最小单元刚度矩阵
k2=BilinearQuadElementStiffness(E2,NU,t,0,0,b,0,b,b,0,b,p);%最小单元刚度矩阵
ke=BilinearQuadElementStiffness(Ee,NU,t,0,0,a,0,a,a,0,a,p);%刚度均匀化等效刚度矩阵
kr=BilinearQuadElementStiffness(Er,NU,t,0,0,a,0,a,a,0,a,p);%柔度均匀化等效刚度矩阵
m1=bilineM(b,b,density1);%质量矩阵
m2=bilineM(b,b,density2);%质量矩阵
Me=bilineM(a,a,density1);%等效后的大质量矩阵
%有限单元法程序段
tnod=zeros(line+1,ar+1);    %有限单元全部节点号存储矩阵
num=line*ar;                %有限单元总数目
%分配节点号程序段
t=0;
for j=1:ar+1    %按列循环
    for i=1:line+1  %按行循环
    t=t+1;
    tnod(i,j)=t;    %为相应单元节点赋值
    end
end 
nodnum=(line+1)*(ar+1); %总节点个数
aixs=zeros(nodnum,2);   %有限元各节点坐标矩阵，第一列为X坐标，第二列为Y坐标
for i=1:ar+1
    aixs(tnod(1:line+1,i),1)=b*(i-1);   %赋x坐标值
end
for i=1:line+1
    aixs(tnod(i,1:ar+1),2)=b*(line+1-i);   %赋y坐标值
end
ZK=billinezk(ar,line,nodnum,tnod,k1,k2);    %组装有限元刚度矩阵
ZM=billinezk(ar,line,nodnum,tnod,m1,m2);    %组装有限元质量矩阵
f=zeros(2*nodnum,1);        %载荷向量,沿X正方向的均布力，总值为 fj
for i=(line+1)*ar+2:nodnum-1    %所加分布载荷的中间节点上分摊到的力
f(2*i-1,1)=fj/(nod-1);
end
f(2*((line+1)*ar+1)-1,1)=fj/(nod-1)/2;  %角点力
f(2*nodnum-1,1)=fj/(nod-1)/2;           %角点力
%**************加边界条件后的有限元总刚、总质量矩阵
ZKC=ZK(2*fstart-1:2*nodnum,2*fstart-1:2*nodnum);    %加边界条件的总刚度矩阵
ZMC=ZM(2*fstart-1:2*nodnum,2*fstart-1:2*nodnum);    %加边界条件的总质量矩阵
%----------有限元位移、频率解-------
u=ZKC\f(2*fstart-1:2*nodnum,1);   %有限元位移解
d(2*fstart-1:2*nodnum,1)=u;
U=0.5*d'*ZK*d;
W=0.5*f'*d;
[v fre]=eig(ZKC,ZMC);   %广义特征值法求解有限元圆频率
ff=zeros(2*freenod,1);  %频率向量
for i=1:2*freenod
ff(i,1)=(fre(i,i)^(1/2))/(2*pi);    %有限元频率解
end
%-----------------绘制有限元位移图、模态图-----
Ux=zeros(nodnum,1);     %全部节点X方向位移
Uy=zeros(nodnum,1);     %全部节点Y方向位移
for i=1:freenod
    Ux(line+1+i,1)=u(2*i-1,1);
    Uy(line+1+i,1)=u(2*i,1);
end
aixsmx=zeros(line+1,ar+1);  %有限元全部节点的X方向坐标矩阵
aixsmy=zeros(line+1,ar+1);  %有限元全部节点的Y方向坐标矩阵
aixsmz=zeros(line+1,ar+1);  %有限元全部节点的Z方向坐标矩阵
aixsmc=zeros(line+1,ar+1);  %网格颜色矩阵
t=1;
for j=1:ar+1
    for i=1:line+1
        aixsmx(i,j)=aixs(t,1)+Ux(t,1);  %原X方向坐标加X方向位移，得到变形后的节点X方向坐标
        aixsmy(i,j)=aixs(t,2)+Uy(t,1);  %原Y方向坐标加Y方向位移，得到变形后的节点Y方向坐标
        aixsmc(i,j)=Uy(t,1);            %全部节点都以Y方向位移为颜色分量
        t=t+1;
    end
end
figure;
hidden off;
mesh(aixsmx,aixsmy,aixsmz,aixsmc);
view(2);
hold on;
%求一个单胞的总刚度矩阵、应变能、外力功
dsubtnod=zeros(dsubline+1,dsubar+1);    %定义单胞内部节点号存储矩阵
%分配单胞内部节点号子程序---------------*******************-----------------
t=0;
for j=1:dsubar+1    %按列循环
    for i=1:dsubline+1  %按行循环
    t=t+1;
    dsubtnod(i,j)=t;    %按由上至下，由左向右的方法，排出特征单元全部节点号
    end
end
dsubnodnum=dsubtnod(dsubline+1,dsubar+1);   %单胞内部总节点个数
%组装一个单胞
zk=billinezk(dsubar,dsubline,dsubnodnum,dsubtnod,k1,k2);
%----有限元法求得的内部胞体位移向量--
dn=zeros(2*dsubnodnum,1);  %内部胞体的变形量
t=1;
for j=1:dsubar+1
    for i=1:dsubline+1
        dn(t,1)=d(2*tnod(dsubline+i,j)-1,1);
        dn(t+1,1)=d(2*tnod(dsubline+i,j),1);
        t=t+2;
    end
end
t_xy=dn(2*(dsubline+1)*dsubar+1,1)/dn(2*(dsubline+1)*dsubar+2,1);     %内部胞体x、y方向位移比例******——————*******
d1=zeros(2*dsubnodnum,1);  %上部胞体的变形量
t=1;
for j=1:dsubar+1
    for i=1:dsubline+1
        d1(t,1)=d(2*tnod(i,j)-1,1);
        d1(t+1,1)=d(2*tnod(i,j),1);
        t=t+2;
    end
end
d3=zeros(2*dsubnodnum,1);  %下部胞体的变形量
t=1;
for j=1:dsubar+1
    for i=1:dsubline+1
        d3(t,1)=d(2*tnod(6+i,j)-1,1);
        d3(t+1,1)=d(2*tnod(6+i,j),1);
        t=t+2;
    end
end
U_ie1=0.5*d1'*zk*d1;    %上部胞体的应变能
U_ie3=0.5*d3'*zk*d3;    %下部胞体的应变能
U_ie=0.5*dn'*zk*dn;     %内部胞体的应变能
Uu=U_ie1+U_ie+U_ie3;    %三个胞体的总应变能
fn=zeros(2*dsubnodnum,1);  %内部胞体所承受的载荷向量
for i=(dsubline+1)*dsubar+2:dsubnodnum-1    %所加分布载荷的中间节点上分摊到的力
fn(2*i-1,1)=fj/subline/dsubline;
end
fn(2*((dsubline+1)*dsubar+1)-1,1)=fj/subline/dsubline/2;            %角点力
fn(2*dsubnodnum-1,1)=fj/subline/dsubline/2;                    %角点力
W_fs=0.5*fn'*dn;
W_fs1=0.5*fn'*d1;
W_fs3=0.5*fn'*d3;
Ww=W_fs1+W_fs+W_fs3;
%一个单元情况下应变能
dc=zk(9:32,9:32)\fn(9:32,1);
d11=zeros(32,1);
d11(9:32)=dc;
U_11=0.5*d11'*zk*d11;
W_11=0.5*fn'*d11;
f1t=zeros(8,1);
f1t(3,1)=fj/6;
f1t(5,1)=f1t(3,1);
dxx=6*U_11/fj;
dyy=dxx/t_xy;
d_xy=zeros(8,1);    %依据能量等效准则，应该产生的位移列向量。
d_xy(3,1)=dxx;
d_xy(4,1)=-dyy;
d_xy(5,1)=dxx;
d_xy(6,1)=dyy;
U_xy=0.5*d_xy'*ke*d_xy;
nxz=U_11/U_xy      %依据能量等效准则，求得的对应于刚度均匀化法求得的刚度矩阵的等效系数。