关于非线性振动IHB法求解过程中的相关问题

wojiuxihuan

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相关程序）：
[size=15.5556px]clear all
clc
tic
%=====输入基本参数（已知条件）===================================
syms t
m=[1.0 1.63 1.0];
c_in=[0.001 0.001 0.001];
Stiffness=[1.0 1.0 1.0];
Stiffness_k=[0.001 0.001 0.001];
F=[0 0 0.1];
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;
%-------------------自由振动--------------------------------------
[vector,value]=eig(K,Mass);     %--------广义特征向量和特征值
[VECTOR,VALUE]=cdf2rdf(vector,value);
omega0=sqrt(VALUE);         
for ii=1:MainNumber
    omega(ii,1)=omega0(ii,ii);%rad/s
end
%===============================================================
Cs=[cos(t) cos(2*t) cos(3*t) sin(t) sin(2*t) sin(3*t)];%见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=[0.1 0.1 0.1 0.1 0.1 0.1]';%%？？？？=====
    A2=[0.1 0.1 0.1 0.1 0.1 0.1]';%%？？？？=====
    A3=[0.1 0.1 0.1 0.1 0.1 0.1]';
for j=1:length(w0)

    A0=[A1;A2;A3];
    %========================三次非线性刚度矩阵====================================
    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=[A01(1),A01(2),A01(3),A01(4),A01(5),A01(6)]';
    A2=[A01(7),A01(8),A01(9),A01(10),A01(11),A01(12)]';
    A3=[A01(13),A01(14),A01(15),A01(16),A01(17),A01(18)]';
end





请教相关问题：1）为什么曲线在共振点处没有出现倾斜？一般在共振处附近一个w对应二个解（稳定解和非稳定解）；2）运用IHB法如何求解非稳定解和稳定解？3）程序在求解过程中是否有错误？

wangsheng37

你算出来都收敛吗？一般非线性幅频曲线都不会是这种，你这看起来像是线性的

Pparis

感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

wojiuxihuan

wangsheng37 发表于 2016-8-3 16:57
你算出来都收敛吗？一般非线性幅频曲线都不会是这种，你这看起来像是线性的

是存在问题，所以发帖向大家咨询！求指导。。。。。。

wojiuxihuan

Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

发出来大家一起看看，希望可以提供帮助，谢谢

Pseudo-lover

楼主问题还没解决吗

wangsheng37

wojiuxihuan 发表于 2016-8-11 09:16
是存在问题，所以发帖向大家咨询！求指导。。。。。。

你的方程是什么？会不会是非线性项太小导致的

Pseudo-lover

wangsheng37 发表于 2016-8-11 14:03
你的方程是什么？会不会是非线性项太小导致的

非线性项过小会导致什么问题

shenyongjun

仔细看，是有倾斜的。响应曲线没有问题，局部放大后效果更好

Edinburgh

shenyongjun 发表于 2016-8-11 21:25
仔细看，是有倾斜的。响应曲线没有问题，局部放大后效果更好

那楼主的问题在哪里？

wojiuxihuan

Pseudo-lover 发表于 2016-8-11 09:51
楼主问题还没解决吗

还没有解决，这是个文献上面的实例（其采用的是增量传递矩阵法）计算结果共振点一致，只是倾斜相差很大。

wojiuxihuan

shenyongjun 发表于 2016-8-11 21:25
仔细看，是有倾斜的。响应曲线没有问题，局部放大后效果更好

不是的，这是个文献上面的实例（其采用的是增量传递矩阵法）计算结果共振点一致，只是倾斜相差很大。

wojiuxihuan

Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

希望大家能提供一些帮助。。。。

wojiuxihuan

Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

希望大家能提供一些帮助。

zghitboy

程序没有报错？