MATLAB有限元仿真出图案例

clear all;
clc;
% for d=1:1:100
syms x y L l1 l2 l3 l4 l5 l6 l7 l8 l9 l10 b h L
% ========材料常数，来源于刘文光轴向共振控制=====
% ====铝合金材料结构尺寸矩形截面%来源于王涛的硕士论文、刘文光材料304不锈钢===陈岩的含裂纹非线性静动力特性，龙慧基于有限元
% L=0.3;%长度单位m
% b=0.02;%宽度
% h=0.02;
A=b*h;%截面积A
I=b*h^3/12;%截面惯性矩
A0=0;%初始损伤的面积
% d=0.05;%直径
% A=0.25*pi*d^2;%截面面积
E=7.3e10;%弹性模量单位Pa
sigmab=5.45e8;%抗拉极限强度
sigmas=4e8;%抗拉强度
rou=2.785e3;%密度
f=668;
mu=0.33;
w=2*pi*f;%圆频率
t=0;
% ======损伤的常数值====
s=2;%一般取值为2-3
SS=125.3e6;%能量损耗值
RR=-1;%对称循环的情况下
DC=0.0001;
% ======振动仿真====
% ====型函数===
LZ=[l1 l2 l3 l4 l5 l6 l7 l8 l9 l10];
% l=L/10;%单元长度分为十个单元
e=x./LZ;
for m=1:10
Nu1=1-e（m）;
Nu2=1-10*e（m）^3+15*e（m）^4-6*e（m）^5;
Nu3=LZ（m）*（e（m）-6*e（m）^3+8*e（m）^4-3*e（m）^5）;
Nu4=LZ（m）^2*（e（m）^2-3*e（m）^3+3*e（m）^4-e（m）^5）/2;
Nu5=e（m）;
Nu6=10*e（m）^3-15*e（m）^4+6*e（m）^5;
Nu7=LZ（m）*（-4*e（m）^3+7*e（m）^4-3*e（m）^5）;
Nu8=LZ（m）^2*（e（m）^3-2*e（m）^4+e（m）^5）/2;
N（::m）=[Nu1 0 0 0 Nu5 0 0 0;
    0 Nu2 Nu3 Nu4 0 Nu6 Nu7 Nu8];
% ------单元质量积分--------
mjf（::m）=eval（int（rou*b*h*N（::m）.‘*N（::m）x0LZ（m）））;%被积分的单元质量矩阵各个单元的质量
% --------------------------
% ====应变转换矩阵，应变与节点位移之间关系采用刚度的定义求解单元刚度矩阵====
Su1=diff（Nu1x）;
Su2=-y*diff（Nu2x2）;
Su3=-y*diff（Nu3x2）;
Su4=-y*diff（Nu4x2）;
Su5=diff（Nu5x）;
Su6=-y*diff（Nu6x2）;
Su7=-y*diff（Nu7x2）;
Su8=-y*diff（Nu8x2）;
S（::m）=[Su1 Su2 Su3 Su4 Su5 Su6 Su7 Su8];
% ------单元刚度矩阵积分-------
 kjf（::m）=E*b*int（int（S（::m）.‘*S（::m）x0LZ（m））y-h/2h/2）;%杆质量矩阵
end

% % % =========动态特性随损伤的变化的情况=====
jj=0;
% D=0;%初始损伤为0
% t=0;
% iii=1;
UUu（:1）=zeros（111）;
% while 1

% % kjf=E*（A-A0）*S.‘*S;
% kjf=E*A*S.‘*S;
% % kjf1=vpa（E*（A-A0）*S.‘*S）;
% % ke=eval（int（kjfx0l））;
% for n=1:10;
%     k（::n）=vpa（eval（int（kjfx0l）））;%杆质量矩阵
% end
% %  k（::1）=vpa（eval（int（kjf1x0l）））;%杆质量矩阵
% ====坐标转换矩阵====
R=[1 0 0 0 0 0 0 0;0 1 0 0 0 0 0 0;0 0 1 0 0 0 0 0;0 0 0 1 0 0 0 0;
    0 0 0 0 1 0 0 0 ;0 0 0 0 0 1 0 0;0 0 0 0 0 0 1 0;0 0 0 0 0 0 0 1 ];
% ====坐标协调矩阵十个单元=====
B1=zeros（844）;
B1（11）=1;B1（22）=1;B1（33）=1;B1（44）=1;B1（55）=1;B1（66）=1;B1（77）=1;B1（88）=1;
B2=zeros（844）;
B2（15）=1;B2（26）=1;B2（37）=1;B2（48）=1;B2（59）=1;B2（610）=1;B2（711）=1;B2（812）=1;
B3=zeros（844）;
B3（19）=1;B3（210）=1;B3（311）=1;B3（412）=1;B3（513）=1;B3（614）=1;B3（715）=1;B3（816）=1;
B4=zeros（844）;
B4（113）=1;B4（214）=1;B4（315）=1;B4（416）=1;B4（517）=1;B4（618）=1;B4（719）=1;B4（820）=1;
B5=zeros（844）;
B5（117）=1;B5（218）=1;B5（319）=1;B5（420）=1;B5（521）=1;B5（622）=1;B5（723）=1;B5（824）=1;
B6=zeros（844）;
B5（121）=1;B5（222）=1;B5（323）=1;B5（424）=1;B6（525）=1;B6（626）=1;B6（727）=1;B6（828）=1;
B7=zeros（844）;
B7（125）=1;B7（226）=1;B7（327）=1;B7（428）=1;B7（529）=1;B7（630）=1;B7（731）=1;B7（832）=1;
B8=zeros（844）;
B8（129）=1;B8（230）=1;B8（331）=1;B8（432）=1;B8（533）=1;B8（634）=1;B8（735）=1;B8（836）=1;
B9=zeros（844）;
B9（133）=1;B9（234）=1;B9（335）=1;B9（436）=1;B9（537）=1;B9（638）=1;B9（739）=1;B9（840）=1;
B10=zeros（844）;
B10（137）=1;B10（238）=1;B10（339）=1;B10（4