关于非线性振动IHB法求解过程中的相关问题
1）研究对象和方法：同样采用上面的三质量非线性系统数据；运用增量谐波平衡法（IHB）法对图1所示的系统进行研究。2）IHB法A1=a1coswt+a2cos2wt+a3cos3wt+b1sinwt+b2sin2wt+b3sin3wt 第一步采用IHB法增量过程；第二部采用IHB法的谐波平衡过程进行求解。3）编写MATLAB程序，经计算可以得到不同w下的a1、a2和a3以及b1、b2和b3（每个w所对应的a1、a2和a3以及b1、b2和b3是唯一的）求解出来的曲线图如图2所示
附（MATLAB相关程序）：
clear all
clc
tic
%=====输入基本参数（已知条件）===================================
syms t
m=;
c_in=;
Stiffness=;
Stiffness_k=;
F=;
Mass=diag(m);%质量矩阵
%=======================内阻尼矩阵=======================================
c=c_in;
MainNumber=3;
%—————————————————————————————————c1
c1(1)=2*c(1); %——1
for ii=2:MainNumber-1
c1(ii)=c(ii)+c(ii-1); %——2 to MainNumber-1
end
c1(MainNumber)=c(MainNumber-1); %——MainNumber
%—————————————————————————————————c2
for ii=1:MainNumber-1
c2(ii)=c(ii);
end
%—————————————————————————————————内阻尼矩阵C
C=diag(c1,0)-diag(c2,1)-diag(c2,-1);
%===============================================================
%========================线性刚度矩阵====================================
k=Stiffness;
MainNumber=3;
%—————————————————————————————————k1
k1(1)=2*k(1); %——1
for ii=2:MainNumber-1
k1(ii)=k(ii)+k(ii-1); %——2 to MainNumber-1
end
k1(MainNumber)=k(MainNumber-1); %——MainNumber
%—————————————————————————————————k2
for ii=1:MainNumber-1
k2(ii)=k(ii);
end
%—————————————————————————————————刚度矩阵K
K=diag(k1,0)-diag(k2,1)-diag(k2,-1);
%K=K*10^6;
%-------------------自由振动--------------------------------------
=eig(K,Mass); %--------广义特征向量和特征值
=cdf2rdf(vector,value);
omega0=sqrt(VALUE);
for ii=1:MainNumber
omega(ii,1)=omega0(ii,ii);%rad/s
end
%===============================================================
Cs=;%见P176公式6.4.8和6.4.9
for i=1:3
for j=1:6
S(i,j+6*(i-1))=Cs(j);
end
end
S1=diff(S,t,1);%对S矩阵进行求1阶导数
S2=diff(S,t,2);%对S矩阵进行求2阶导数
w0=0.01:0.1:2;
%==================假定初始值============================================
A1=';%%？？？？=====
A2=';%%？？？？=====
A3=';
for j=1:length(w0)
A0=;
%========================三次非线性刚度矩阵====================================
q0=S*A0;
q00=q0*q0';
No_linear_k=Stiffness_k*q00;
nonlinear_k=No_linear_k;
MainNumber=3;
%—nonlinear————————————————————————————————k1
nonlinear_k1(1)=2*nonlinear_k(1); %——1
for ii=2:MainNumber-1
nonlinear_k1(ii)=nonlinear_k(ii)+nonlinear_k(ii-1); %——2 to MainNumber-1
end
nonlinear_k1(MainNumber)=nonlinear_k(MainNumber-1); %——MainNumber
%—————————————————————————————————k2
for ii=1:MainNumber-1
nonlinear_k2(ii)=nonlinear_k(ii);
end
%—————————————————————————————————刚度矩阵K
nonlinear_K=diag(nonlinear_k1,0)-diag(nonlinear_k2,1)-diag(nonlinear_k2,-1);
%K=K*10^6;
%===============================================================
%======================M=========================================
fm=inline(S'*Mass*S2);
MM=quadv(fm,0,2*pi);%见P176公式6.4.15，M
%======================C=========================================
fc=inline(S'*C*S1);
CC=quadv(fc,0,2*pi);%见P176公式6.4.15，M
%=====================K==========================================
fk=inline(S'*K*S);
KK=quadv(fk,0,2*pi);
%=====================nonlinear_K==========================================
fk_nonlinear=inline(S'*nonlinear_K*S);
nonlinear_KK=quadv(fk_nonlinear,0,2*pi);
%===============================================================
fF=inline(S'*F'*sin(t));
FF=quadv(fF,0,2*pi);
%=============w0=1以及A0=========================================
Kmc=w0(j)^2*MM+w0(j)*CC+KK+3*nonlinear_KK;
R=FF-(w0(j)^2*MM+w0(j)*CC+KK+nonlinear_KK)*A0;
Rmc=-(w0(j)^2*MM+CC)*A0;
%===============================================================
deta_A=inv(Kmc)*(R+Rmc*0); %drtA1(1)
A01=A0+deta_A;
n=1;
tol=1;
epsR=0.01;
while tol>epsR
A0=A01;
%======================M=========================================
fm=inline(S'*Mass*S2);
MM=quadv(fm,0,2*pi);%见P176公式6.4.15，M
%======================C=========================================
fc=inline(S'*C*S1);
CC=quadv(fc,0,2*pi);%见P176公式6.4.15，M
%=====================K==========================================
fk=inline(S'*K*S);
KK=quadv(fk,0,2*pi);
%=====================nonlinear_K==========================================
fk_nonlinear=inline(S'*nonlinear_K*S);
nonlinear_KK=quadv(fk_nonlinear,0,2*pi);
%===============================================================
fF=inline(S'*F'*sin(t));
FF=quadv(fF,0,2*pi);
%===============================================================
%%%%%带入推导出的公式
Kmc=w0(j)^2*MM+w0(j)*CC+KK+3*nonlinear_KK;
R=FF-(w0(j)^2*MM+w0(j)*CC+KK+nonlinear_KK)*A0;
Rmc=-(w0(j)^2*MM+CC)*A0;
%%%%%
tol=norm(R);
if(n>1000)
disp('迭代步数太多，可能不收敛')
return;
end
%===============================================================
deta_A=inv(Kmc)*(R+Rmc*0); %drtA1(1)
A01=A0+deta_A;
n=n+1;
end
for i=1:length(A01)
AA(j,i)=A01(i);
end
A1=';
A2=';
A3=';
end
请教相关问题：1）为什么曲线在共振点处没有出现倾斜？一般在共振处附近一个w对应二个解（稳定解和非稳定解）；2）运用IHB法如何求解非稳定解和稳定解？3）程序在求解过程中是否有错误？
你算出来都收敛吗？一般非线性幅频曲线都不会是这种，你这看起来像是线性的 感谢楼主分享经验程序写的很清晰可以看出楼主做研究很严谨 值得学习 wangsheng37 发表于 2016-8-3 16:57
你算出来都收敛吗？一般非线性幅频曲线都不会是这种，你这看起来像是线性的
是存在问题，所以发帖向大家咨询！求指导。。。。。。 Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验程序写的很清晰可以看出楼主做研究很严谨 值得学习
发出来大家一起看看，希望可以提供帮助，谢谢 楼主问题还没解决吗 wojiuxihuan 发表于 2016-8-11 09:16
是存在问题，所以发帖向大家咨询！求指导。。。。。。
你的方程是什么？会不会是非线性项太小导致的 wangsheng37 发表于 2016-8-11 14:03
你的方程是什么？会不会是非线性项太小导致的
非线性项过小会导致什么问题 仔细看，是有倾斜的。响应曲线没有问题，局部放大后效果更好 shenyongjun 发表于 2016-8-11 21:25
仔细看，是有倾斜的。响应曲线没有问题，局部放大后效果更好
那楼主的问题在哪里？ Pseudo-lover 发表于 2016-8-11 09:51
楼主问题还没解决吗
还没有解决，这是个文献上面的实例（其采用的是增量传递矩阵法）计算结果共振点一致，只是倾斜相差很大。 shenyongjun 发表于 2016-8-11 21:25
仔细看，是有倾斜的。响应曲线没有问题，局部放大后效果更好
不是的，这是个文献上面的实例（其采用的是增量传递矩阵法）计算结果共振点一致，只是倾斜相差很大。 Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验程序写的很清晰可以看出楼主做研究很严谨 值得学习
希望大家能提供一些帮助。。。。 Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验程序写的很清晰可以看出楼主做研究很严谨 值得学习
希望大家能提供一些帮助。 程序没有报错？