[已答复] 增量谐波平衡法程序

[复制链接]
阿灿啊 发表于 2019-2-25 11:19:35
根据网上那个范德波极限环程序编的,为什么我的幅频图是很奇怪的折线


clear all
close all
clc
tic
%=====输入基本参数(已知条件)===================================
%t=1;
syms t
syms w0
m0=0.2712e6;k0=0.3707e4;c0=0.4439e4;v21=1.6039;
mu1=0.05;mu2=0.005;
m1=m0*mu1;m2=m0*mu2;
xi1=0.1195;xi2=0.2177;%阻尼比
lambda1=0.9266;lambda2=0.8117;%频率比
w00=sqrt(k0/m0);
w1=lambda1*w00;w2=lambda2*w00;
k1=w1^2*m1;
k2=w1^2*m2;
c1=2*xi1*m1*w1;
c2=2*xi2*m2*w2;
%=======================质量矩阵m=======================================
m=diag([m0/v21 m1 m2]);
%=======================阻尼矩阵c=======================================
c=[c0/v21+c1*v21+c2*v21 -c1*v21 -c2*v21;
    -c1*v21 c1*v21 0;
    -c2*v21 0 c2*v21];
%========================线性刚度矩阵k====================================
k=[k0/v21+k1*v21 -k1*v21 0;
    -k1*v21 k1*v21 0;
    0 0 0];
%========================外激励f====================================
f=[0.5836;0;0]*1e6;
%===============================================================
Cs=[cos(t) cos(3*t) cos(5*t) sin(t) sin(3*t) sin(5*t)];
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阶导数
%==================假定初始值============================================
A1=[1 1 1 1 1 1]';
A2=[1 1 1 1 1 1]';
A3=[1 1 1 1 1 1]';
A0=[A1;A2;A3];
%========================三次非线性刚度矩阵kn====================================
nonlinear_k=[k2*v21 0 -k2*v21;0 0 0;-k2*v21 0 k2*v21];
x0=S*A0;
q00=(x0(3)-x0(1))^2;
kn=nonlinear_k*q00;
%======================M=========================================
fm=inline(S'*m*S2);
M=quadv(fm,0,2*pi);
%======================C=========================================
fc=inline(S'*c*S1);
C=quadv(fc,0,2*pi);
%=====================K==========================================
fk=inline(S'*k*S);
K=quadv(fk,0,2*pi);
%=====================nonlinear_K==========================================
fk_nonlinear=inline(S'*kn*S);
KN=quadv(fk_nonlinear,0,2*pi);
%=============================F==================================
fF=inline(S'*f*cos(t));
F=quadv(fF,0,2*pi);
w0=0.4;
%=============w0=1以及A0=========================================
Kmc=w0^2*M+w0*C+K+3*KN;
R=F-(w0^2*M+w0*C+K+KN)*A0;
Rmc=-(w0^2*M+C)*A0;
%===============================================================
deta_A=inv(Kmc)*(R+Rmc*0);     %drtA1(1)
A01=A0+deta_A;

n=1;
tol=1;
epsR=0.001;

while tol>epsR
     A0=A01;
%======================M=========================================
fm=inline(S'*m*S2);
M=quadv(fm,0,2*pi);
%======================C=========================================
fc=inline(S'*c*S1);
C=quadv(fc,0,2*pi);
%=====================K==========================================
fk=inline(S'*k*S);
K=quadv(fk,0,2*pi);
%=====================nonlinear_K==========================================
fk_nonlinear=inline(S'*kn*S);
KN=quadv(fk_nonlinear,0,2*pi);
%=============================F==================================
fF=inline(S'*f*cos(t));
F=quadv(fF,0,2*pi);
    %===============================================================
    %%%%%带入推导出的公式
Kmc=w0^2*M+w0*C+K+3*KN;
R=F-(w0^2*M+w0*C+K+KN)*A0;
Rmc=-(w0^2*M+C)*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

     X0=S*A0;
     dX0=S1*A0;
     figure(1)
     plot(A01);hold on;grid on;
tt=0:.1:100;
x01=subs(X0(1),tt);
x02=subs(X0(2),tt);
x03=subs(X0(3),tt);
dx01=subs(dX0(1),tt);
dx02=subs(dX0(2),tt);
dx03=subs(dX0(3),tt);
figure(2)
plot(x01,dx01,'b','linewidth',2)
hold on
plot(x02,dx02,'b','linewidth',2)
hold on
plot(x03,dx03,'b','linewidth',2)
xlabel('x0')

4 条回复


阿灿啊 发表于 2019-2-25 11:42:21
本帖最后由 阿灿啊 于 2019-2-26 21:40 编辑

帖子有误,删不了。。。

untitled.jpg

wx_hXfN5DY4 发表于 2019-12-6 09:16:19
请问您这个问题解决了吗?我想请教一下关于谐波平衡法的程序

Yu_P 发表于 2020-12-10 22:24:59
你这个不是幅频曲线

博博士 发表于 2021-4-8 23:23:56
这个没有幅频曲线啊,某bao搜索“增量谐波平衡法”,还有弧长法以及各种IHB的推广
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
热门教程
站长推荐
快速回复 返回顶部 返回列表