p080123xz3.asv

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

%p080123xz3.m  验证精度，对不同胞体进行计算。
%将位移、频率结果进行比较，完备的程序，验证正确！
%——————基本变量——————
clear;clc;
E1=200e3;E2=E1;     %第一种情况：均匀材料，各向同性。
%E1=200e3;E2=E1/50; %第二种情况：最角上硬，中间软
%E2=200e3;E1=E2/50; %第三种情况：最角上软，中间硬
t=1;NU=0.3;p=1;      %弹性模量、厚度、泊松比、平面问题的种类，1为平面应力，2为应变
density1=5.0e-6;density2=5.0e-6;fj=100;
%————————————————————等效系数——————————
tt=0.976234073823371;    %由p080123.m文件求出等效系数，nxy，拷贝入本程序，赋值给tt.
%----------------****************------------------**********-------------
line=15;ar=30;  %原始有限元的单元行数，列数
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=5;subar=10;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);   %各节点坐标矩阵
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);    %载荷向量,上下两侧同时受均布力，方向相反
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=zeros(32,1);
d(2*fstart-1:2*nodnum,1)=u;
[v fre]=eig(ZKC,ZMC);   %有限元圆频率解
ff=zeros(2*freenod,1);  % (nodnum-line-1),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);
aixsmy=zeros(line+1,ar+1);
aixsmz=zeros(line+1,ar+1);
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);
        aixsmy(i,j)=aixs(t,2)+Uy(t,1);
        aixsmc(i,j)=Uy(t,1);
        t=t+1;
    end
end
%-----------------------------------
%特征单元法
dsubtnum=zeros(8,1); %定义当前组装等效单元节点在总刚中的对应行号存储程序
dsubttnum=zeros(2,2); %定义当前组装等效单元节点号存储程序
dsubtnod=zeros(dsubline+1,dsubar+1);%定义总节点号存储矩阵
dsubnum=dsubline*dsubar;%单胞中单元个数
%分配等效单元节点号子程序---------------*******************-----------------
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);   %总节点个数
kg=zeros(8,8);
kg1=zeros(8,8);
kg1=ke; 
kg=ke*tt;
%组装特征单元
subtnum=zeros(8,1); %定义当前组装等效单元节点在总刚中的对应行号存储程序
subttnum=zeros(2,2); %定义当前组装等效单元节点号存储程序
subtnod=zeros(subline+1,subar+1);%定义总节点号存储矩阵
subnum=subline*subar;%总单胞个数
t=0;
for j=1:subar+1    %按列循环
    for i=1:subline+1  %按行循环
    t=t+1;
    subtnod(i,j)=t;%排出单元节点号
    end
end
subnodnum=subtnod(subline+1,subar+1);   %总节点个数
%--------设定等效节点坐标-------------**********************-------------
subaixs=zeros(subnodnum,2);   %各节点坐标矩阵
for i=1:subar+1
    subaixs(subtnod(1:subline+1,i),1)=(i-1)*a;   %赋x坐标值
end
for i=1:subline+1
    subaixs(subtnod(i,1:subar+1),2)=(subline+1-i)*a;   %赋y坐标值
end
%-----------组装等效后单元的刚度、质量矩阵------------**************--------------
ZKg=zeros(subnodnum*2,subnodnum*2);    %组装全部单元刚度矩阵
ZKe=zeros(subnodnum*2,subnodnum*2);    %组装全部单元刚度矩阵
ZKr=zeros(subnodnum*2,subnodnum*2);    %组装全部单元刚度矩阵
ZMe=zeros(subnodnum*2,subnodnum*2);    %组装全部等效质量矩阵
for j=1:subar
    for i=1:subline
    subttnum=subtnod(i:i+1,j:j+1);
    subtnum(1,1)=2*subttnum(2,1)-1;
    subtnum(2,1)=2*subttnum(2,1);
    subtnum(3,1)=2*subttnum(2,2)-1;
    subtnum(4,1)=2*subttnum(2,2);
    subtnum(5,1)=2*subttnum(1,2)-1;
    subtnum(6,1)=2*subttnum(1,2);
    subtnum(7,1)=2*subttnum(1,1)-1;
    subtnum(8,1)=2*subttnum(1,1);
    ZKg=Asube(ZKg,kg,subtnum); %特征等效法总等效刚度矩阵
    ZKe=Asube(ZKe,ke,subtnum); %刚度等效法等效总刚度矩阵
    ZKr=Asube(ZKr,kr,subtnum); %柔度等效法等效总刚度矩阵    
    ZMe=Asube(ZMe,Me,subtnum); %总等效质量矩阵    
    end
end
fg=zeros(2*subnodnum,1);    %等效载荷
for i=(subar*(subline+1)+2):subnodnum-1
    fg(2*i-1,1)=fj/(subnod-1);
end
fg(2*((subline+1)*subar+1)-1,1)=fj/(subnod-1)/2;
fg(2*subnodnum-1,1)=fj/(subnod-1)/2;
ug=ZKg(2*(subline+2)-1:2*subnodnum,2*(subline+2)-1:2*subnodnum)\fg(2*(subline+2)-1:2*subnodnum,1);  %特征等效位移解
ue=ZKe(2*(subline+2)-1:2*subnodnum,2*(subline+2)-1:2*subnodnum)\fg(2*(subline+2)-1:2*subnodnum,1);  %刚度等效位移解
ur=ZKr(2*(subline+2)-1:2*subnodnum,2*(subline+2)-1:2*subnodnum)\fg(2*(subline+2)-1:2*subnodnum,1);  %柔度度等效位移解
de=zeros(2*subnodnum,1);
de(2*subfstart-1:2*subnodnum,1)=ue;
U_ave=0.5*de'*ZKe*de;
ng=a/max(max(ug))/3;
ne=a/max(max(ue))/3;
nr=a/max(max(ur))/3;
%------------------绘制等效位移图------------************------------
uxe=zeros(subnodnum,1);
uye=zeros(subnodnum,1);
uxg=zeros(subnodnum,1);
uyg=zeros(subnodnum,1);
uxr=zeros(subnodnum,1);
uyr=zeros(subnodnum,1);
for i=1:subfreenod
    uxe(subline+1+i,1)=ue(2*i-1,1);
    uye(subline+1+i,1)=ue(2*i,1);
    uxg(subline+1+i,1)=ug(2*i-1,1);
    uyg(subline+1+i,1)=ug(2*i,1);
    uxr(subline+1+i,1)=ur(2*i-1,1);
    uyr(subline+1+i,1)=ur(2*i,1);
end
%----变形网格------
subaixsmxe=zeros(subline+1,subar+1);
subaixsmye=zeros(subline+1,subar+1);
subaixsmze=zeros(subline+1,subar+1);
subaixsmce=zeros(subline+1,subar+1);
subaixsmxg=zeros(subline+1,subar+1);
subaixsmyg=zeros(subline+1,subar+1);
subaixsmzg=zeros(subline+1,subar+1);
subaixsmcg=zeros(subline+1,subar+1);
subaixsmxr=zeros(subline+1,subar+1);
subaixsmyr=zeros(subline+1,subar+1);
subaixsmzr=zeros(subline+1,subar+1);
subaixsmcr=zeros(subline+1,subar+1);
t=1;
for j=1:subar+1
    for i=1:subline+1
        subaixsmxe(i,j)=subaixs(t,1)+uxe(t,1);
        subaixsmye(i,j)=subaixs(t,2)+uye(t,1);
        subaixsmce(i,j)=uxe(t,1);
        subaixsmxg(i,j)=subaixs(t,1)+uxg(t,1);
        subaixsmyg(i,j)=subaixs(t,2)+uyg(t,1);
        subaixsmcg(i,j)=uxg(t,1);
        subaixsmxr(i,j)=subaixs(t,1)+uxr(t,1);
        subaixsmyr(i,j)=subaixs(t,2)+uyr(t,1);
        subaixsmcr(i,j)=uxr(t,1);
        t=t+1;
    end
end
%-------求解模态、频率
[vE freE]=eig(ZKe(2*subfstart-1:2*subnodnum,2*subfstart-1:2*subnodnum),ZMe(2*subfstart-1:2*subnodnum,2*subfstart-1:2*subnodnum));   %有限元圆频率解
[vG freG]=eig(ZKg(2*subfstart-1:2*subnodnum,2*subfstart-1:2*subnodnum),ZMe(2*subfstart-1:2*subnodnum,2*subfstart-1:2*subnodnum));   %有限元圆频率解
[vR freR]=eig(ZKr(2*subfstart-1:2*subnodnum,2*subfstart-1:2*subnodnum),ZMe(2*subfstart-1:2*subnodnum,2*subfstart-1:2*subnodnum));   %有限元圆频率解
tmax=max(max(vE));
Tmax=max(max(v));
ffe=zeros(2*subfreenod,1);
ffg=zeros(2*subfreenod,1);
ffr=zeros(2*subfreenod,1);
for i=1:2*subfreenod  % i 的取值范围取决于加没有被约束的自由度数目
    ffe(i,1)=(freE(i,i)^(1/2))/(2*pi);    %有限元频率解
    ffg(i,1)=(freG(i,i)^(1/2))/(2*pi);    %有限元频率解
    ffr(i,1)=(freR(i,i)^(1/2))/(2*pi);    %有限元频率解
end
[ffe1,se1]=sort(ffe);
[ffg1,sg1]=sort(ffg);
[ffr1,sr1]=sort(ffr);
err1=zeros(20,1);
err2=zeros(20,1);
err3=zeros(20,1);
for i=1:20
    err1(i,1)=abs((ffe1(i,1)-ff(i,1))/ff(i,1));
    err2(i,1)=abs((ffg1(i,1)-ff(i,1))/ff(i,1));
    err3(i,1)=abs((ffr1(i,1)-ff(i,1))/ff(i,1));
end
figure;
plot(err1,'s','color','w');
hold on;
plot(err2,'.','color','w');
plot(err3,'o','color','w');
legend('UNIS-50','Present-50','UNIC-50');
xlabel('No. of Frequency');
ylabel('Relative of Error');
merr1=max(err1);
merr2=max(err2);
figure;
plot(ff(1:20,1),'+','color','w');
hold on
plot(ffe1(1:20,1),'s','color','w');
plot(ffg1(1:20,1),'.','color','w');
plot(ffr1(1:20,1),'o','color','w');
legend('FEM-450','UNIS-50','Present-50','UNIC-50');
xlabel('No. of Frequency');
ylabel('Frequency');
hold off