对我如下模型进行Campbell图的绘制,模型使用转轴为timoshenko梁单元,以及建立好了mm、kk、cc、gg矩阵。 close all;clc;clear;
%% 输入轴/输出轴结构参数: [phi直径 length长度 分段数]
Es=2.07e11;%弹性模量 Pa(N/m^2)
Poissons=0.29;
Gs=Es/(2*(1+Poissons));%剪切模量
dens=7840; %材料密度 kg/m3
%%输入轴
para1=[58*2*1e-3 60*2*1e-3 2;
58*2*1e-3 40*1e-3 1;
58*2*1e-3 20*2*1e-3 2];
%%输出轴
para2=[100*2*1e-3 30*2*1e-3 2;
100*2*1e-3 15*2*1e-3 2];
%%
%%单元信息 :[单元编号,节点1编号,节点2编号,直径,弹性模量,剪切模量]
num=0;
%%输入轴
for i=1:length(para1(:,1))
seg=para1(i,3); %分段数
for j=1:seg
num=num+1;
ele(num,:)=[num num num+1];
end
end
%%输出轴
for i=1:length(para2(:,1))
seg=para2(i,3); %分段数
Ds=para2(i,1); %直径
for j=1:seg
num=num+1;
ele(num,:)=[num num+1 num+2 ];
end
end
%%弹簧支承单元 单元信息[单元号 节点1 节点2 kx ky kz] ele_num 一根梁划分单元数
num = num+1;
ele(num,:) = [num 2 num+2 ];%第一个节点弹簧支座单元
num = num+1;
ele(num,:) = [num 3 num+2 ];%第一根梁第三个节点
num = num+1;
ele(num,:) = [num 7 num+2 ];%第二根梁 第7个节点
num = num+1;
ele(num,:) = [num 11 num+2 ];%第二根梁 第11个节点
%%
%%节点位置:[节点编号,x,y,z]
num=0;
%%输入轴
for i=1:length(para1(:,1))
seg=para1(i,3);
len=para1(i,2);
if i==1
len_befor=0;
else
len_befor=sum(para1(1:i-1,2));
end
for j=1:seg
if i==1 && j==1
num=num+1;
nodeCoord(num,:)=[num,0,0,0];%节点位置:[节点编号,x,y,z]
num=num+1;
nodeCoord(num,:)=[num,0,0,len_befor+len/seg*j];
else
num=num+1;
nodeCoord(num,:)=[num,0,0,len_befor+len/seg*j];
end
end
end
%%输出轴
for i=1:length(para2(:,1))
seg=para2(i,3); %分段数
len=para2(i,2); %分段长度
if i==1
len_befor=0;
else
len_befor=sum(para2(1:i-1,2));
end
for j=1:seg
if i==1 && j==1
num=num+1;
nodeCoord(num,:)=[num,-0.09,0,0.2];%节点位置:[节点编号,x,y,z]
num=num+1;
nodeCoord(num,:)=[num,-0.09+len_befor+len/seg*j,0,0.2];
else
num=num+1;
nodeCoord(num,:)=[num,-0.09+len_befor+len/seg*j,0,0.2];
end
end
end
% 支承弹簧节点信息[节点位置:[节点编号,x,y,z]
num = num + 1;
nodeCoord(num,:)=[num,nodeCoord(1,2),nodeCoord(1,3)-100,nodeCoord(1,4)];%第一根梁第1个节点
num = num + 1;
nodeCoord(num,:)=[num,nodeCoord(3,2),nodeCoord(3,3)-100,nodeCoord(3,4)];%第一根梁倒数第3个节点
num = num + 1;
nodeCoord(num,:)=[num,nodeCoord(7,2),nodeCoord(7,3)-100,nodeCoord(7,4)];%第二根梁 第7个节点
num = num + 1;
nodeCoord(num,:)=[num,nodeCoord(11,2),nodeCoord(11,3)-100,nodeCoord(11,4)];%第二根梁 第11个节点
%%
node_num=size(nodeCoord,1); %节点数
ele_num=size(ele,1); %轴段单元数
BarLength=BarsLength(nodeCoord(:,2)',nodeCoord(:,3)',nodeCoord(:,4)',ele); %计算每个单元的长度
xx = nodeCoord(:,2);
yy = nodeCoord(:,3);
zz = nodeCoord(:,4);
%% 齿轮参数设置
Ep = 2.07e11; %弹性模量
vp = 0.29;%泊松比
Gp = Ep*0.5/(1+vp);%剪切模量
denp = 7840;
Z1 = 23;Z2 = 143;
%小齿轮半径pinion 齿数23
Rp = 0.002*25/2*cosd(25);
Rp_in = 0.058;
Rg = 0.28;
Rg_in = 0.245;
Wp = 0.04;%齿宽
Wg = 0.035;
%% 齿轮动态参数设置*********************************
% 扭矩
Tg = 950;
Tp = Tg*Z1/Z2;
%%%%%***********转速
N1 =5000;%r/min 小齿轮转速
N2 = N1*Z2/Z1;
%
fr = N1/60;%主动轮转频
wp = 2*pi*fr;%主动轮转角速度rad/s ,匀速
wg = wp*Z2/Z1;%从动轮转角速度rad/s,匀速
omegap = wp*Z1;%主动轮啮合角速度rad/s
%% beam K M G
beam_K = zeros(12, 12, ele_num); % 创建用于存储 beam_K 的三维数组
beam_M = zeros(12, 12, ele_num);
beam_G = zeros(12, 12, ele_num);
%%输入轴
Num_e1=sum(para1(:,3));
Ds = para1(1,1); % 输入轴直径
for i = 1:Num_e1
[beam_Ki, beam_Mi, beam_Gi] = get_beam(Es, Gs, BarLength(i),Poissons, Ds, 0.05*2, 10/9, dens);%通过子函数get_beam_K计算轴段单元的M/C/K矩阵。
%[beam_K, beam_M, beam_G] = get_beam_K(弹性模量, 剪切模量, 单元长度, 单元外径, 单元内径, 校正因子, 材料密度)
beam_K(:, :, i) = beam_Ki; % 将每次迭代得到的 beam_Ki 存储在三维数组中的对应位置
beam_M(:, :, i) = beam_Mi;
beam_G(:, :, i) = beam_Gi; %__________________________获得梁单元的质量、刚度矩阵、阻尼矩阵
end
%%输出轴
Ds = para2(1,1); % 输出轴直径
for i = (Num_e1+1):ele_num
[beam_Ki, beam_Mi, beam_Gi] = get_beam(Es, Gs, BarLength(i),Poissons, Ds, 0.09*2, 10/9, dens);%通过子函数get_beam_K计算轴段单元的M/C/K矩阵。
%[beam_K, beam_M, beam_G] = get_beam_K(弹性模量, 剪切模量, 单元长度, 单元外径, 单元内径, 校正因子, 材料密度)
beam_K(:, :, i) = beam_Ki; % 将每次迭代得到的 beam_Ki 存储在三维数组中的对应位置
beam_M(:, :, i) = beam_Mi;
beam_Gi = beam_Gi*wg/wp;
beam_G(:, :, i) = beam_Gi; %__________________________获得梁单元的质量、刚度矩阵、阻尼矩阵
end
%% gear M G
[gear_Mp, gear_Gp] = get_gear(denp,Wp,Rp,0.058*2);
[gear_Mg, gear_Gg] = get_gear(denp,Wg,Rg,0.245*2);
gear_M = [gear_Mp,zeros(6,6);zeros(6,6),gear_Mg];
gear_G = [gear_Gp,zeros(6,6);zeros(6,6),gear_Gg*wg/wp];
%% supportint K C
%四个支承弹簧
% [bearing_K1,bearing_C1] = get_bearing(8e5, 8e5, 7e5, 7e5, 7e5, 0,4e2, 4e2, 4e2, 1e2, 1e2, 0);
% [bearing_K2,bearing_C2] = get_bearing(8e5, 8e5, 7e5, 7e5, 7e5, 0,4e2, 4e2, 4e2, 1e2, 1e2, 0);
% [bearing_K3,bearing_C3] = get_bearing(8e5, 8e5, 7e5, 7e5, 7e5, 0,4e2, 4e2, 4e2, 1e2, 1e2, 0);
% [bearing_K4,bearing_C4] = get_bearing(8e5, 8e5, 7e5, 7e5, 7e5, 0,4e2, 4e2, 4e2, 1e2, 1e2, 0);
%约束:节点编号,x方向约束,y方向约束,z方向约束,theta_x theta_y theta_z
constr=[12 0 0 0 0 0 0;
13 0 0 0 0 0 0;
14 0 0 0 0 0 0;
15 0 0 0 0 0 0];
%% 支承弹簧刚度
%1 3 节点支承弹簧
kx1 = 1.24e8;ky1 = 1.6e8;kz1 = 4.37e6;
kthetax1 = 4.18e3;kthetay1 = 1.06e3;kthetaz1 = 0;
%7 11 节点支承弹簧
kx2 = 1.59e8;ky2 = 2.02e8;kz2 = 3.11e6;
kthetax2 = 7.72e3;kthetay2 = 4.43e3;kthetaz2 = 0;
%% 组装带集中质量圆盘TIMOSHENKO梁 M K G 矩阵
%建立空白矩阵
Dofs = 6*node_num;
Ke=zeros(Dofs,Dofs); %总体刚度矩阵
Me=zeros(Dofs,Dofs); %总体质量矩阵
Ce=zeros(Dofs,Dofs); %总体阻尼矩阵
Ge=zeros(Dofs,Dofs); %总体陀螺矩阵
%遍历所有单元,将各单元阵分块组装到总体阵
for iEle =1:9
if iEle<=5
%该单元的两个节点的编号
n1=ele(iEle,2);n2=ele(iEle,3);
%计算坐标变换矩阵
O = zeros(3,3);
T = [1 0 0;
0 1 0;
0 0 1];
R=[T O O O;
O T O O;
O O T O;
O O O T];
%计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
ke= R'*beam_K(:, :, iEle)*R;
me = R'*beam_M(:, :, iEle)*R;
ge = R'*beam_G(:, :, iEle)*R;
%将各单元刚度分块组装到总刚相应位置
eleDof=[n1*6-5:n1*6,n2*6-5:n2*6];
Ke(eleDof,eleDof)=Ke(eleDof,eleDof)+ke;
Me(eleDof,eleDof)=Me(eleDof,eleDof)+me;
Ge(eleDof,eleDof)=Ge(eleDof,eleDof)+ge;
else
%该单元的两个节点的编号
n1=ele(iEle,2);n2=ele(iEle,3);
%计算坐标变换矩阵
O = zeros(3,3);
T = [1 0 0;
0 0 1;
0 1 0];
R=[T O O O;
O T O O;
O O T O;
O O O T];
%计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
ke= R'*beam_K(:, :, iEle)*R;
me = R'*beam_M(:, :, iEle)*R;
ge = R'*beam_G(:, :, iEle)*R;
%将各单元刚度分块组装到总刚相应位置
eleDof=[n1*6-5:n1*6,n2*6-5:n2*6];
Ke(eleDof,eleDof)=Ke(eleDof,eleDof)+ke;
Me(eleDof,eleDof)=Me(eleDof,eleDof)+me;
Ge(eleDof,eleDof)=Ge(eleDof,eleDof)+ge;
%轴承单元
% %该单元的两个节点的编号
% n1=ele(iEle,2);n2=ele(iEle,3);
% %计算坐标变换矩阵
% T_bearing = [1 0 0;
% 0 0 -1;
% 0 1 0];
% R_bearing=[T_bearing O O O;
% O T_bearing O O;
% O O T_bearing O;
% O O O T_bearing];
% %计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
% ke_bearing= get_bearing(kx1,ky1,kz1,kthetax1,kthetay1,kthetaz1);
% % 将各单元刚度分块组装到总刚相应位置
% ke_bearing = R_bearing'*ke_bearing*R_bearing;
% eleDof=[n1*6-5:n1*6,n2*6-5:n2*6];
% Ke(eleDof,eleDof)=Ke(eleDof,eleDof)+ke_bearing;
end
end
%% 轴系瑞丽阻尼
a00 = 0;
a11 = 1e-7;
C=a00*Me+a11*Ke; %形成阻尼矩阵
Ce = Ce + C;
%%
meshingnode1 = 5;%啮合单元节点
meshingnode2 = 10;
%组装刚性圆盘质量矩阵
%%啮合自由度
meshingdof = [meshingnode1*6-5:meshingnode1*6,meshingnode2*6-5:meshingnode2*6];
% n1=meshingnode1;n2=meshingnode2;
%计算坐标变换矩阵
T_meshing = [1 0 0;
0 0 -1;
0 -1 0];
R_mehsing=[T_meshing O O O;
O T_meshing O O;
O O T_meshing O;
O O O T_meshing];
Me(meshingdof,meshingdof) = Me(meshingdof,meshingdof) + R_mehsing'*gear_M*R_mehsing;
Ge(meshingdof,meshingdof) = Ge(meshingdof,meshingdof) + R_mehsing'*gear_G*R_mehsing;
%% 带集中质量圆盘+弹性支承系统 瑞丽阻尼
%%
%添加约束 在四弹簧支座单元处添加 固定约束
for i = 1:length(constr(:,1))
restraineddof(1+(i-1)*6:6+(i-1)*6) = (constr(i,1)-1)*6+(1:6);
end
activeDof = setdiff((1:Dofs)',restraineddof);
%% 支承轴承
%轴承单元
for iEle =10:13
if iEle<=11
%该单元的两个节点的编号
n1=ele(iEle,2);n2=ele(iEle,3);
%计算坐标变换矩阵
T_bearing = [0 0 1;
1 0 0;
0 -1 0];
R_bearing=[T_bearing O O O;
O T_bearing O O;
O O T_bearing O;
O O O T_bearing];
%计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
ke_bearing= get_bearing(kx1,ky1,kz1,kthetax1,kthetay1,kthetaz1);
ce_bearing1 = [1000 0 0 0 0 0;
0 1000 0 0 0 0;
0 0 1000 0 0 0;
0 0 0 1000 0 0;
0 0 0 0 1000 0;
0 0 0 0 0 0];
ce_bearing = [ce_bearing1 ,-ce_bearing1;
-ce_bearing1, ce_bearing1];
% 将各单元刚度分块组装到总刚相应位置
ke_bearing = R_bearing'*ke_bearing*R_bearing;
ce_bearing = R_bearing'*ce_bearing*R_bearing;
eleDof=[n1*6-5:n1*6,n2*6-5:n2*6];
Ke(eleDof,eleDof)=Ke(eleDof,eleDof)+ke_bearing;
Ce(eleDof,eleDof)=Ce(eleDof,eleDof)+ce_bearing;
else
%该单元的两个节点的编号
n1=ele(iEle,2);n2=ele(iEle,3);
%计算坐标变换矩阵
T_bearing = [0 0 1;
1 0 0;
0 -1 0];
R_bearing=[T_bearing O O O;
O T_bearing O O;
O O T_bearing O;
O O O T_bearing];
%计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
ke_bearing= get_bearing(kx2,ky2,kz2,kthetax2,kthetay2,kthetaz2);
ce_bearing1 = [1000 0 0 0 0 0;
0 1000 0 0 0 0;
0 0 1000 0 0 0;
0 0 0 1000 0 0;
0 0 0 0 1000 0;
0 0 0 0 0 0];
ce_bearing = [ce_bearing1 ,-ce_bearing1;
-ce_bearing1, ce_bearing1];
% 将各单元刚度分块组装到总刚相应位置
ke_bearing = R_bearing'*ke_bearing*R_bearing;
ce_bearing = R_bearing'*ce_bearing*R_bearing;
eleDof=[n1*6-5:n1*6,n2*6-5:n2*6];
Ke(eleDof,eleDof)=Ke(eleDof,eleDof)+ke_bearing;
Ce(eleDof,eleDof)=Ce(eleDof,eleDof)+ce_bearing;
end
end
%% 时变啮合刚度 —— 取平均值近似
kt_mean = 4.1182e8; % 你的固定值
% kt_mean = 0;
% 啮合方向投影向量 W(12×1 列向量)
rp = 39.1e-3; rg = 243e-3;
alpha = 25*pi/180;
fai = 90*pi/180; %#ok<NASGU>
% 你的 V_proj 是行向量,这里改成列向量并外积 —— FIX #2
W = [-cos(alpha); -sin(alpha); 0; 0; 0; rp; ...
cos(alpha); sin(alpha); 0; -rg*sin(alpha); -rg*cos(alpha); 0];
% 啮合刚度矩阵 Km = kt * W*W^T(12×12)
Km_local = kt_mean * (W * W');
% 变换到全局并装配 —— 与文献式(32)一致
Ke(meshingdof, meshingdof) = Ke(meshingdof, meshingdof) + R_mehsing' * Km_local * R_mehsing;
%% 施加“地”约束:把 4 个地节点完全固定(通过删减自由度)
% 这四个地节点就是你构造的地基弹簧下端节点:编号 12~15(按你前面拓扑生成)
constr = [12 0 0 0 0 0 0;
13 0 0 0 0 0 0;
14 0 0 0 0 0 0;
15 0 0 0 0 0 0];
% restraineddof = [];
% for i = 1:size(constr,1)
% restraineddof = [restraineddof, (constr(i,1)-1)*6 + (1:6)]; %#ok<AGROW>
% end
% restraineddof = unique(restraineddof);
% activeDof = setdiff((1:ndof)', restraineddof);
%添加约束 在四弹簧支座单元处添加 固定约束
for i = 1:length(constr(:,1))
restraineddof(1+(i-1)*6:6+(i-1)*6) = (constr(i,1)-1)*6+(1:6);
end
activeDof = setdiff((1:Dofs)',restraineddof);
% —— 关键:形成“已约束”的缩减矩阵 —— FIX #1
mm = Me(activeDof, activeDof);
kk = Ke(activeDof, activeDof);
cc = Ce(activeDof, activeDof);
gg = Ge(activeDof, activeDof);
最新发布