查看: 3204|回复: 36

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

[复制链接]
发表于 2016-7-20 10:02 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
1)研究对象和方法:同样采用上面的三质量非线性系统数据;运用增量谐波平衡法(IHB)法对图1所示的系统进行研究。
2)IHB法
A1=a1coswt+a2cos2wt+a3cos3wt+b1sinwt+b2sin2wt+b3sin3wt
    第一步采用IHB法增量过程;第二部采用IHB法的谐波平衡过程进行求解。
   3)编写MATLAB程序,经计算可以得到不同w下的a1a2a3以及b1b2b3(每个w所对应的a1a2a3以及b1b2b3是唯一的)求解出来的曲线图如图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)程序在求解过程中是否有错误?

图1非线性系统参数图

图1非线性系统参数图

图2 计算结果

图2 计算结果
回复
分享到:

使用道具 举报

发表于 2016-8-3 16:57 | 显示全部楼层
你算出来都收敛吗?一般非线性幅频曲线都不会是这种,你这看起来像是线性的
发表于 2016-8-4 08:49 | 显示全部楼层
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习
 楼主| 发表于 2016-8-11 09:16 | 显示全部楼层
wangsheng37 发表于 2016-8-3 16:57
你算出来都收敛吗?一般非线性幅频曲线都不会是这种,你这看起来像是线性的

是存在问题,所以发帖向大家咨询!求指导。。。。。。
 楼主| 发表于 2016-8-11 09:24 | 显示全部楼层
Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

发出来大家一起看看,希望可以提供帮助,谢谢
发表于 2016-8-11 09:51 | 显示全部楼层
楼主问题还没解决吗
发表于 2016-8-11 14:03 | 显示全部楼层
wojiuxihuan 发表于 2016-8-11 09:16
是存在问题,所以发帖向大家咨询!求指导。。。。。。

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

点评

非线性项过小会导致什么问题  详情 回复 发表于 2016-8-11 14:19
发表于 2016-8-11 14:19 | 显示全部楼层
wangsheng37 发表于 2016-8-11 14:03
你的方程是什么?会不会是非线性项太小导致的

非线性项过小会导致什么问题
发表于 2016-8-11 21:25 | 显示全部楼层
仔细看,是有倾斜的。响应曲线没有问题,局部放大后效果更好

点评

那楼主的问题在哪里?  详情 回复 发表于 2016-8-12 08:51
发表于 2016-8-12 08:51 | 显示全部楼层
shenyongjun 发表于 2016-8-11 21:25
仔细看,是有倾斜的。响应曲线没有问题,局部放大后效果更好

那楼主的问题在哪里?
 楼主| 发表于 2016-8-17 09:46 | 显示全部楼层
Pseudo-lover 发表于 2016-8-11 09:51
楼主问题还没解决吗

还没有解决,这是个文献上面的实例(其采用的是增量传递矩阵法)计算结果共振点一致,只是倾斜相差很大。
 楼主| 发表于 2016-8-17 09:46 | 显示全部楼层
shenyongjun 发表于 2016-8-11 21:25
仔细看,是有倾斜的。响应曲线没有问题,局部放大后效果更好

不是的,这是个文献上面的实例(其采用的是增量传递矩阵法)计算结果共振点一致,只是倾斜相差很大。
 楼主| 发表于 2016-8-17 09:49 | 显示全部楼层
Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

希望大家能提供一些帮助。。。。
 楼主| 发表于 2016-8-17 09:52 | 显示全部楼层
Pparis 发表于 2016-8-4 08:49
感谢楼主分享经验  程序写的很清晰  可以看出楼主做研究很严谨 值得学习

希望大家能提供一些帮助。
发表于 2016-8-18 08:30 | 显示全部楼层
程序没有报错?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表