增量谐波平衡法（IHB）简介，附Matlab程序
增量谐波平衡法在国内的应用非常广泛，已经应用到简单梁、运动梁、框架、板壳、叶片、齿轮系统、轴系系统、时滞系统等很多系统。增量谐波平衡法的兄弟姐妹也很多，但基本原理都是一样的，根据具体应用做了适宜的发展。关键的地方是谐波系数的确定。谐波系数可以公式算，可以用FFT程序，可以数值积分，也可以最小二乘拟合，等等。国内有学者用抗混频的FFT来算，也是很好的改进。还有学者修改了谐波假设，考虑了线性趋势项。现在的算法主要是采用正余弦函数，最后是实数线性方程组。如果采用复数指数函数，得到是复数线性方程组。两个方法显然是等价的。
更复杂的应用是计算频响曲线，稳定性，分岔响应等，可以拾阶而上。另外，虽然增量谐波平衡法主要应用在响应计算，也可以应用到求解逆问题。比方非线性系统的参数识别，已经有低自由度系统的数值验证和实验应用。
——以上来自科学网suguang
方法的推导
一般地，工程结构中非线性振动问题可以归结为如下的方程：
式中：m为质量; c为阻尼;k(q)为切线刚度矩阵;q为位移，f 为荷载。
令：
变量置换：
则：
则式(2)变换为：
令：
代入式(3)中，得：
即：
其中：
用傅里叶级数展开与时间有关的各项：
代入式(4)中并应用谐波平衡法，得到：
其中
而
残余力矢
上式中
到此，得到了一个增量形式的非线性代数方程组(5)，如果选择一个初始点(ω0，q)， 给定一个Δω，通过迭代可得到{Δq}。
——以上摘自熊铁华，常晓林发表于《武汉大学学报(工学版)》的《非线性振动系统的主共振的增量谐波平衡法》一文。
实例：两个耦合范德波振子
具体可参考唐南发表于中山大学研究生专刊自然科学版的《应用于范德波方程的增量谐波平衡法》一文。
Matlab代码：
clear allclctic%定义各参数syms tw0=3;epsR=0.001;m=;epsilon=0.24;r=0.4;delta=0.56;k=;Cs=;Cs1=diff(Cs，t，1);S=;A1=';A2=';A0=;T1=;T2=;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=;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，求wwa1=0.0;%变换矩阵，使ww变量代替a1Kmc11=-Rmc(:，1);Kmcr=;%求未知变量AA=inv(Kmcr)*R;%drtA1(1)ww=AA(1);%drtW(1)%赋予新变量新值A01=A0+;%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>epsRA0=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;endKmc11=-Rmc(:，1);Kmcr=;
AA=inv(Kmcr)*R;ww=AA(1);%A00=;A01=A0+;w01=w0+ww;n=n+1;endX0=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 onplot(xo2，dxo2，'b'，'linewidth'，2)axis([-3 3 -3 3])title('范德波极限环')xlabel('x0')ylabel('dx0')toc
运行结果：范德波极限环
——以上代码由声振之家会员zhangwenjing分享，代码未经验证。
很不错！！！ 增量谐波平衡法matlab程序，单自由度，多自由度均可以，弧长增量法求解多值区间，不稳定解，获得复杂幅频曲线，并判断解的稳定性和分岔。可以参考TB店：https://item.taobao.com/item.htm ... ;abbucket=19#detail 博博士 发表于 2019-12-11 17:41
增量谐波平衡法matlab程序，单自由度，多自由度均可以，弧长增量法求解多值区间，不稳定解，获得复杂幅频曲 ...
增量谐波平衡法matlab程序，单自由度多自由度均可以，弧长增量法求解多值区间、不稳定解，追踪获得复杂幅频曲线，并判断周期解的稳定性和分岔，并可与数值解验证，结果吻合良好。可以参考TB店铺：增量谐波平衡法，shop517677650.taobao.com，QQ973199830 IHBM的确有效，但是如下两个问题足够你吃一壶：
1，初值不好设计；初值选取不当，不收敛，或者自己设计的收敛准则，收敛以后解无意义；
2，对测试数据要求较高，比如你要位移，速度（导数）和加速度测试结果，但是对MDOF系统，你很难做到；
mxlzhenzhu 发表于 2020-1-27 13:04
IHBM的确有效，但是如下两个问题足够你吃一壶：
1，初值不好设计；初值选取不当，不收敛，或者自己设计的 ...
初始选取都是根据经验猜想，也可以用线性解作为初值，或者利用其它方法先近似算一下，然后设置与精确解相近的初始，收敛就容易多了。测试数据要求较高？对于多自由度系统，无非是计算时间更长一些而已，尤其是在想要得到精度较高解的时候，需要采用较多的谐波项数。
页:
[1]