clear all
clc
tic
%=====输入基本参数(已知条件)===================================
%t=1;
syms t
syms w0
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];
% % % m1=[0.0213,1.3307,0.0509,0.3727];
% % % Cs=[cos(t) cos(3*t) sin(t) sin(3*t)];%见P176公式6.4.8和6.4.9
% % % Cs1=diff(Cs,t,1);%对Cs矩阵进行求1阶导数
% % % KK=[Cs,Cs];
% % % SS=diag(KK);%=[cos(t) cos(3*t) sin(t) sin(3*t)];%见P176公式6.4.8和6.4.9
% % % 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)];%见P176公式6.4.10
% % % S1=S';
Mass=diag(m);%质量矩阵或转动惯量矩阵
%=======================内阻尼矩阵=======================================
c=c_in;
MainNumber=3;
%—————————————————————————————————c1
c1(1)=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)=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;
%===============================================================
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阶导数
%==================假定初始值============================================
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];
%========================三次非线性刚度矩阵====================================
q0=S*A0;
q00=q0*q0';
No_linear_k=Stiffness_k*q00;
nonlinear_k=No_linear_k;
MainNumber=3;
%—nonlinear————————————————————————————————k1
nonlinear_k1(1)=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=0.4;
%=============w0=1以及A0=========================================
Kmc=w0^2*MM+w0*CC+KK+3*nonlinear_KK;
R=FF-(w0^2*MM+w0*CC+KK+nonlinear_KK)*A0;
Rmc=-(w0^2*MM+CC)*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'*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^2*MM+w0*CC+KK+3*nonlinear_KK;
R=FF-(w0^2*MM+w0*CC+KK+nonlinear_KK)*A0;
Rmc=-(w0^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
X0=S*A0;
%dX0=S1*A0;
没有合适的资源?快使用搜索试试~ 我知道了~
基于matlab实现的增量谐波平衡方法(IHB)是一种有效的定量分析方法
共1个文件
m:1个
1.该资源内容由用户上传,如若侵权请联系客服进行举报
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
2.虚拟产品一经售出概不退款(资源遇到问题,请及时私信上传者)
版权申诉
0 下载量 159 浏览量
2024-05-02
15:23:53
上传
评论
收藏 2KB RAR 举报
温馨提示
基于matlab实现的增量谐波平衡方法(IHB)是一种有效的定量分析方法,既适用于分析弱非线性振动系统,也适用于分析强非线性振动系统,并且可以提高周期解的精度.rar
资源推荐
资源详情
资源评论
收起资源包目录
基于matlab实现的增量谐波平衡方法(IHB)是一种有效的定量分析方法,既适用于分析弱非线性振动系统,也适用于分析强非线性振动系统,并且可以提高周期解的精度.rar (1个子文件)
IHB.m 6KB
共 1 条
- 1
资源评论
依然风yrlf
- 粉丝: 891
- 资源: 3118
下载权益
C知道特权
VIP文章
课程特权
开通VIP
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功