增量谐波平衡法（IHB）简介，附Matlab程序

weixin

　　增量谐波平衡法在国内的应用非常广泛，已经应用到简单梁、运动梁、框架、板壳、叶片、齿轮系统、轴系系统、时滞系统等很多系统。增量谐波平衡法的兄弟姐妹也很多，但基本原理都是一样的，根据具体应用做了适宜的发展。

　　关键的地方是谐波系数的确定。谐波系数可以公式算，可以用FFT程序，可以数值积分，也可以最小二乘拟合，等等。国内有学者用抗混频的FFT来算，也是很好的改进。还有学者修改了谐波假设，考虑了线性趋势项。现在的算法主要是采用正余弦函数，最后是实数线性方程组。如果采用复数指数函数，得到是复数线性方程组。两个方法显然是等价的。

　　更复杂的应用是计算频响曲线，稳定性，分岔响应等，可以拾阶而上。另外，虽然增量谐波平衡法主要应用在响应计算，也可以应用到求解逆问题。比方非线性系统的参数识别，已经有低自由度系统的数值验证和实验应用。

　　——以上来自科学网suguang

　　方法的推导
　　一般地，工程结构中非线性振动问题可以归结为如下的方程：

　　式中：m为质量; c为阻尼;k(q)为切线刚度矩阵;q为位移，f 为荷载。
　　令：

　　变量置换：

　　则：

　　则式(2)变换为：

　　令：

　　代入式(3)中，得：

　　即：

　　其中：

　　用傅里叶级数展开与时间有关的各项：

　　代入式(4)中并应用谐波平衡法，得到：

　　其中

　　而

　　残余力矢

　　上式中

　　到此，得到了一个增量形式的非线性代数方程组(5)，如果选择一个初始点(ω0，q)， 给定一个Δω，通过迭代可得到{Δq}。

　　——以上摘自熊铁华，常晓林发表于《武汉大学学报(工学版)》的《非线性振动系统的主共振的增量谐波平衡法》一文。

　　实例：两个耦合范德波振子
　　具体可参考唐南发表于中山大学研究生专刊自然科学版的《应用于范德波方程的增量谐波平衡法》一文。

　　Matlab代码：

clear all

clc

tic

%定义各参数

syms t

w0=3;

epsR=0.001;

m=[1 0；0 1];

epsilon=0.24;

r=0.4;

delta=0.56;

k=[1+epsilon*r -epsilon*r;-epsilon*r 1+epsilon*delta];

Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];

Cs1=diff(Cs，t，1);

S=[cos(t) cos(3*t) sin(t) sin(3*t) 0 0 0 0;0 0 0 0 cos(t) cos(3*t) sin(t) sin(3*t)];

A1=[1 1 1 1]';

A2=[1 1 1 1]';

A0=[A1；A2];

T1=[eye(4，4) zeros(4，4)];

T2=[zeros(4，4) eye(4，4)];

S2=diff(S，t，2);

fm=inline(S'*m*S2);

M=quadv(fm，0，2*pi);

fk=inline(S'*k*S);

K=quadv(fk，0，2*pi);

S1=diff(S，t，1);

c=[1 0；0 1];

fc=inline(S'*c*S1);

C=quadv(fc，0，2*pi);

c3=diag(S*A0).^2;

fc3=inline(S'*c3*S1);

C3=quadv(fc3，0，2*pi);

k2=diag(S*A0).*diag(S1*A0);

fk2=inline(S'*k2*S);

K2=quadv(fk2，0，2*pi);

%代入推导出的公式

Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;

R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;

Rmc=-(2*w0*M+epsilon*(C3-C))*A0;

%AA首元素已知a1=0.0，求ww

a1=0.0;

%变换矩阵，使ww变量代替a1

Kmc11=-Rmc(:，1);

Kmcr=[Kmc11 Kmc(:，2:size(Kmc，2))];

%求未知变量

AA=inv(Kmcr)*R;

%drtA1(1)

ww=AA(1);

%drtW(1)

%赋予新变量新值

A01=A0+[a1; AA(2:length(A0)，1)];

%A(1)+drtA(1)

% Aw0=AA+A00;

%A1(0)+drtA1(1)=A1(1)

w01=w0+ww;

%W+drtW(1)

n=1;

tol=1;

while tol>epsR

A0=A01;

w0=w01;

c3=diag(S*A0).^2;

fc3=inline(S'*c3*S1);

C3=quadv(fc3，0，2*pi);

k2=diag(S*A0).*diag(S1*A0);

fk2=inline(S'*k2*S);

K2=quadv(fk2，0，2*pi);

%带入推导出的公式

Kmc=w0^2*M+epsilon*w0*(C3-C)+K+2*epsilon*w0*K2;

R=-(w0^2*M+epsilon*w0*(C3-C)+K)*A0;

Rmc=-(2*w0*M+epsilon*(C3-C))*A0;

%%%%%

tol=norm(R);

if(n>1000)

disp('迭代步数太多，可能不收敛')

return;

end

Kmc11=-Rmc(:，1);

Kmcr=[Kmc11 Kmc(:，2:size(Kmc，2))];

AA=inv(Kmcr)*R;

ww=AA(1);

%A00=[w0；A0(2:6，1)];

A01=A0+[a1；AA(2:length(A0)，1)];

w01=w0+ww;

n=n+1;

end

X0=S*A0;

dX0=S1*A0;

%绘范德波图

tt=0:.1:10;

xo1=subs(X0(1)，tt);

xo2=subs(X0(2)，tt);

dxo1=subs(dX0(1)，tt);

dxo2=subs(dX0(2)，tt);

figure(1)

plot(xo1，dxo1，'b'，'linewidth'，2)

hold on

plot(xo2，dxo2，'b'，'linewidth'，2)

axis([-3 3 -3 3])

title('范德波极限环')

xlabel('x0')

ylabel('dx0')

toc

　　运行结果：范德波极限环

　　——以上代码由声振之家会员zhangwenjing分享，代码未经验证。

zhiyuan2008

很不错！！！

博博士

增量谐波平衡法matlab程序，单自由度，多自由度均可以，弧长增量法求解多值区间，不稳定解，获得复杂幅频曲线，并判断解的稳定性和分岔。可以参考TB店：https://item.taobao.com/item.htm ... ;abbucket=19#detail

博博士

博博士 发表于 2019-12-11 17:41
增量谐波平衡法matlab程序，单自由度，多自由度均可以，弧长增量法求解多值区间，不稳定解，获得复杂幅频曲 ...

增量谐波平衡法matlab程序，单自由度多自由度均可以，弧长增量法求解多值区间、不稳定解，追踪获得复杂幅频曲线，并判断周期解的稳定性和分岔，并可与数值解验证，结果吻合良好。可以参考TB店铺：增量谐波平衡法，shop517677650.taobao.com，QQ973199830

mxlzhenzhu

IHBM的确有效，但是如下两个问题足够你吃一壶：
1，初值不好设计；初值选取不当，不收敛，或者自己设计的收敛准则，收敛以后解无意义；
2，对测试数据要求较高，比如你要位移，速度（导数）和加速度测试结果，但是对MDOF系统，你很难做到；

花花三公子

mxlzhenzhu 发表于 2020-1-27 13:04
IHBM的确有效，但是如下两个问题足够你吃一壶：
1，初值不好设计；初值选取不当，不收敛，或者自己设计的 ...

初始选取都是根据经验猜想，也可以用线性解作为初值，或者利用其它方法先近似算一下，然后设置与精确解相近的初始，收敛就容易多了。测试数据要求较高？对于多自由度系统，无非是计算时间更长一些而已，尤其是在想要得到精度较高解的时候，需要采用较多的谐波项数。