请教用Newmark法求振动方程的响应的问题
下面是用ode求解振动响应的m文件，由于例的方程是有限元振动方程，想把下面的程序用Newmark法求解响应只用一个m来求解，不知道怎么弄，请高手指点！主程序
:
function BallBrg_NonL_Forum
% 求解外圈固定球轴承的变柔度(VC-Varying Compliance)振动（基于赵凌燕的论文）
clear
clc
%% 参数设置
% 用了全局变量来传递一些变量，不推荐，但是懒得改了，好心人优化一下。
global w d D Nb gama kn M C F
% 为了方便绘制分岔图而设置的参数
n_One_T = 100;% 每个周期的采样点数
n_T = 100;% 采样时间占几个周期
% 61903/P5(17*30*7) 球轴承参数
d=0.0173;% 内滚道直径
D=0.0265;% 外滚道直径
Nb=9; % 滚子数
n_n = 0;
w_limit1=100;% 最低转速（rpm）
w_limit2=20000;% 最高转速（rpm）
w_step = 100;% 转速变化步长(rpm)
q_initial(1:4,1) = 1e-11;% 初始值
gama = 0.00002;% 间隙(m)
F = 6;% 径向力(N)
kn = 7.055e9*0.001^1.5;
% 滚子与滚道之间接触力与变形量的关系(N/mm^1.5)。赵的论文给出。
M=0.6*;% 质量矩阵
C=200*;% 阻尼矩阵
%% 响应计算循环
for w_rpm=w_limit1:w_step:w_limit2
n_n = n_n+1 % 计数变量
disp(w_rpm)
w = w_rpm*pi/30;% 转化为rad/s单位
wi = w;% 内圈角速度
wo = 0;% 外圈角速度
w_cage = ( wi*d/2+wo*D/2 )/2/((D+d)/4);% 保持架
w_vc = w_cage*Nb/2/pi; % 变刚度频率(vc频率)。单位Hz
T_vc = 1/w_vc;% vc周期
dt=T_vc/n_One_T;% 取点时间步长，dt随转速变化。
time=n_T*T_vc;% 总的时间
n = round(time/dt);% 离散点数
t_span(1:n) = linspace(0,time,n);% 时间数组
= ode23tb('BallBrg_NonL_Sub_Forum', t_span, q_initial);
% 至于用什么ode函数求解合适需要比较验证
Q{n_n}=q;
save Q.mat Q; % 存储数据
end
disp('Calculation is done!')
子程序
[ - ]CODE:
function dq = BallBrg_NonL_Sub_Forum(t,q)
% BallBrg_NonL调用的微分方程子程序
global w d D Nb gama kn M C F
wi = w;
wo = 0;
w_cage=( wi*d+wo*D )/4/((D+d)/4);% 保持架转速（rad/s）
fq=zeros(2,1);% 轴承力初值
diff_1_3 = q(1,1);% 水平方向位移
diff_2_4 = q(2,1);% 垂直方向位移
% 求轴承的非线性反力
for No_ball=1:Nb
sita(No_ball) = 2*pi/Nb*(No_ball-1) + w_cage*t;% 第No_ball个滚珠的位置角
Clearance(No_ball,1) = diff_1_3*sin( sita(No_ball) ) ...
+ diff_2_4*cos( sita(No_ball) ) - gama;% 滚珠与内滚道的间隙变化。
% 判断哪几个滚动体受到接触力
if Clearance(No_ball)<=0;
Clearance(No_ball) = 0;
end
fs = abs( (1000*Clearance(No_ball))^1.5 );
fq(1,1) = fq(1,1)+kn*fs*sin(sita(No_ball));
fq(2,1) = fq(2,1)+kn*fs*cos(sita(No_ball));
end
F_m1d1_cos = 0;% 不平衡力在水平方向的投影。本例不考虑。
F_m1d1_sin = 0;% 不平衡力在垂直方向的投影。本例不考虑。
Fq(1,1)= - fq(1,1) + F_m1d1_cos;% 水平方向外力
Fq(2,1)= - fq(2,1) + F_m1d1_sin - F;% 垂直方向外力
K = ;% 刚性转子，轴段为刚性。
% 动力学微分方程
dq(3:4,1)=inv(M)*(Fq-K*q(1:2,1)-C*q(3:4,1));% x和y方向加速度
dq(1:2,1)=q(3:4,1); 由于是用有限元列的振动方程，质量阵，刚度阵比较大用ｏｄｅ４５求解比较麻烦，在网上看到有用Newmark求响应的，用ｏｄｅ４５求响应时需要用的主程序和子程序两个m文件，用Newmark求时是不是只要一个m函数文件就可以了 ，程序中的循环比较多不知道怎么求解
请大家讨论～～ 下面是徐荣桥老师书中exam8_2例子中的一段，不知对你有用否
step3. 计算时程响应(Newmark法)
% step3.1 初始计算
gama = 0.5 ;
beta = 0.25 ;
C = zeros( size( gK ) ) ;
= size( gK ) ;
alpha0 = 1/beta/gDeltaT^2 ;
alpha1 = gama/beta/gDeltaT ;
alpha2 = 1/beta/gDeltaT ;
alpha3 = 1/2/beta - 1 ;
alpha4 = gama/beta - 1 ;
alpha5 = gDeltaT/2*(gama/beta-2) ;
alpha6 = gDeltaT*(1-gama) ;
alpha7 = gama*gDeltaT ;
K1 = gK + alpha0*gM + alpha1*C;
timestep = floor(gTimeEnd/gDeltaT) ;
% step3.2 对K1进行边界条件处理
= size( gBC1 ) ;
K1im = zeros(N,bc1_number) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
K1im(:,ibc) = K1(:,m) ;
K1(:,m) = zeros( node_number*3, 1 ) ;
K1(m,:) = zeros( 1, node_number*3 ) ;
K1(m,m) = 1.0 ;
end
= lu(K1) ; % 进行三角分解，节省后面的求解时间
% step3.3 计算初始加速度
gAcce(:,1) = gM\(-gK*gDisp(:,1)-C*gVelo(:,1)) ;
% step3.4 对每一个时间步计算
for i=2:1:timestep
if mod(i,100) == 0
fprintf( '当前时间步：%d\n', i ) ;
end
f1 = gM*(alpha0*gDisp(:,i-1)+alpha2*gVelo(:,i-1)+alpha3*gAcce(:,i-1)) ...
+ C*(alpha1*gDisp(:,i-1)+alpha4*gVelo(:,i-1)+alpha5*gAcce(:,i-1)) ;
% 对f1进行边界条件处理
= size( gBC1 ) ;
for ibc=1:1:bc1_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*3 + d ;
f1 = f1 - gBC1(ibc,3) * K1im(:,ibc) ;
f1(m) = gBC1(ibc,3) ;
end
y = KL\f1 ;
gDisp(:,i) = KU\y ;
gAcce(:,i) = alpha0*(gDisp(:,i)-gDisp(:,i-1)) - alpha2*gVelo(:,i-1) - alpha3*gAcce(:,i-1) ;
gVelo(:,i) = gVelo(:,i-1) + alpha6*gAcce(:,i-1) + alpha7*gAcce(:,i) ;
end 谢谢楼主！我看看
回复 #3 shanlihong 的帖子
能不能说下这本书的具体名字？我想详细看下。多谢 结构分析的有限元与MATLAB程序设计 收藏了，回去好好研究一下。谢谢楼主。
页:
[1]