Demo entry 6749297

111

   

Submitted by 111 on Jun 12, 2018 at 16:07
Language: Matlab. Code size: 2.9 kB.

clc;clear;
l1=0.05;l2=0.065;l3=0.05;l4=0.045;
r1=0.01;r2=0.075;r3=0.06;r4=0.05;
p=7850;   
E=2090.00;%单位MPa 强度系数
A=0.25*pi*E;
for n=1:4   %产生四个m
    eval(['m',num2str(n),'=pi*p*(r',num2str(n),'^2)',';']);
end

for n=1:4   %产生四个Jz
    eval(['J',num2str(n),'=m',num2str(n),'*r',num2str(n),'^2',';']);
end

% % % % % % % %产生四个R% % % % % % % % % 
R=zeros(4,10);   
for n=1:4   
    eval(['R',num2str(n),'=R',';']);
end
R1=eye(4,10);
R2(:,3:6) =R1(:,1:4);
R3(:,5:8) =R2(:,3:6);
R4(:,7:10) =R3(:,5:8);
% % % % % % % %把数据带入集中质量矩阵 % % % % % % % % % 
M=zeros(4,4);
for n=1:4   %产生四个M,4x4零矩阵
    eval(['M',num2str(n),'=M',';']);
end

for n=1:4   %把数据带入集中质量矩阵M
    
    eval(['M',num2str(n),'(1,1)=m',num2str(n),'*l',num2str(n),';']);
    eval(['M',num2str(n),'(2,2)=J',num2str(n),';']);
    eval(['M',num2str(n),'(3,3)=m',num2str(n),'*l',num2str(n),';']);
    eval(['M',num2str(n),'(4,4)=J',num2str(n),';']);
end

for n=1:4   %产生10x10质量矩阵取名MMi
    eval(['MM',num2str(n),'=R',num2str(n),'''*M',num2str(n),'*R',num2str(n),';']);
end
MM=MM1+MM2+MM3+MM4;            %将4个集中质量矩阵MMi求和得到总质量阵MM

% % % % % % % % % % % %计算K刚度矩阵 % % % % % % % % % % % 
K=zeros(4,4);
for n=1:4   %产生四个K,4x4零矩阵
    eval(['K',num2str(n),'=K',';']);
end

for n=1:4   
    %对角线4个元素先写好
    eval(['K',num2str(n),'(1,1)=12;']);
    eval(['K',num2str(n),'(2,2)=4*l',num2str(n),'^2',';']);
    eval(['K',num2str(n),'(3,3)=12',';']);
    eval(['K',num2str(n),'(4,4)=4*l',num2str(n),'^2',';']);
    %写右上角元素
    eval(['K',num2str(n),'(1,2)=6*l',num2str(n),';']);
    eval(['K',num2str(n),'(1,3)=-12;']);
    eval(['K',num2str(n),'(1,4)=6*l',num2str(n),';']);
    eval(['K',num2str(n),'(2,3)=-6*l',num2str(n),';']);
    eval(['K',num2str(n),'(2,4)=2*l',num2str(n),'^2',';']);
    eval(['K',num2str(n),'(3,4)=-6*l',num2str(n),';']);
end
     %将上三角矩阵对称一下,得到4个大K
K1 = triu(K1,0) + tril(K1',-1); 
K2 = triu(K2,0) + tril(K2',-1);
K3 = triu(K3,0) + tril(K3',-1);
K4 = triu(K4,0) + tril(K4',-1);


for n=1:4   %四个K乘以系数EI/l^3
    eval(['K',num2str(n),'=A*(r',num2str(n),'^4)/(l',num2str(n),'^3)*K',num2str(n),';']);
end

for n=1:4   %产生10x10刚度矩阵取名KKi
    eval(['KK',num2str(n),'=R',num2str(n),'''*K',num2str(n),'*R',num2str(n),';']);
end
KK=KK1+KK2+KK3+KK4;            %将4个集中质量矩阵KKi求和得到总刚度阵KK
% % % % % % % % % % % % % % % % % % % % % % % % % % % % 

t = 100;
F = 200 * sin(t) + 150 * cos(2 * t);

Q = [ 0 0 0 0 0 0 0 0 0 F]';
dq = @(t,q)[q(9:16); MM \ (-KK * q(1:8) + Q)]; %动力学方程
q0 = zeros(16,1);
[t,y] = ode45(dq, [0 1], q0);
plot(t,y(:,1:8)); %位移坐标变化曲线

N = 4;
[V, D] =eig((MM) \ KK); %求解特征向量
for i = 1 : N;
    omg(i) = sqrt(D(i, i)); %求解前四阶固有频率
end;
t = (0 : 0.1 : 2);
x = (0 : 0.0125 : 0.25);
[t, x] = meshgrid(t, x);
y = 32 * x * (200 * sin(t) + 150 * cos(2 * t)) / (pi * 0.01^3);
mesh(x, t, y); %悬臂梁受力曲线图
xlabel('x')
ylabel('t')
zlabel('应力')

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).