| 查看: 1702 | 回复: 0 | ||
[求助]
用Newmark方法计算系统的动力学响应的matlab程序
|
|
请大家帮忙看看这个程序有什么问题？用Newmark方法计算系统的动力学响应，结果大的惊人。 function[Q,V,AA]=newmarkb E=2.1e11;P=7850;D1=0.405;d1=0.375;D2=0.375;d2=0.335;D3=0.335;d3=0.285;D4=0.285;d4=0.225;D5=0.225;d5=0.150; A=(pi*(D1^2-d1^2))/4; I=(pi*(D1^4-d1^4))/64; M1= Mass (P,A,I,0,0,13,0); A=(pi*(D2^2-d2^2))/4; I=(pi*(D2^4-d2^4))/64; M2= Mass (P,A,I,13,0,26,0); A=(pi*(D3^2-d3^2))/4; I=(pi*(D3^4-d3^4))/64; M3= Mass (P,A,I,26,0,39,0); A=(pi*(D4^2-d4^2))/4; I=(pi*(D4^4-d4^4))/64; M4= Mass (P,A,I,39,0,52,0); A=(pi*(D5^2-d5^2))/4; I=(pi*(D5^4-d5^4))/64; M5= Mass (P,A,I,52,0,65,0); M=zeros(18,18); M= MAssemble(M,M1,1,2); M= MAssemble(M,M2,2,3); M= MAssemble(M,M3,3,4); M= MAssemble(M,M4,4,5); M= MAssemble(M,M5,5,6);% 整体质量矩阵 A=(pi*(D1^2-d1^2))/4; I=(pi*(D1^4-d1^4))/64; K1= LStiffness1 (E,A,I,0,0,13,0); A=(pi*(D2^2-d2^2))/4; I=(pi*(D2^4-d2^4))/64; K2= LStiffness2 (E,A,I,13,0,26,0); A=(pi*(D3^2-d3^2))/4; I=(pi*(D3^4-d3^4))/64; K3= LStiffness3 (E,A,I,26,0,39,0); A=(pi*(D4^2-d4^2))/4; I=(pi*(D4^4-d4^4))/64; K4= LStiffness4 (E,A,I,39,0,52,0); A=(pi*(D5^2-d5^2))/4; I=(pi*(D5^4-d5^4))/64; K5= LStiffness5 (E,A,I,52,0,65,0); K=zeros(18,18); K= KAssemble(K,K1,1,2); K= KAssemble(K,K2,2,3); K= KAssemble(K,K3,3,4); K= KAssemble(K,K4,4,5); K= KAssemble(K,K5,5,6);% 整体刚度矩阵 C=0.00776*M+0.00398*K; %瑞利阻尼矩阵 %newmark系数 dt=0.01; nt=20;%计算相应时间 betae=0.25; alfa=0.5; a0=1/(betae*dt^2); a1=alfa/(betae*dt); a2=1/(betae*dt); a3=1/(2*betae)-1; a4=alfa/betae-1; a5=dt/2*(alfa/betae-2); a6=dt*(1-alfa); a7=dt*alfa; %系数定义完 KE=K+a0*M+a1*C; [L,U]=lu(KE); q=zeros(18,1); v=zeros(18,1); pp=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-100,0,100].'; a=M^(-1)*(pp-K*q-C*v); t=0; Q(:,1)=q; V(:,1)=v; AA(:,1)=a; PP(:,1)=pp; for i=1:nt-1 PP(:,i+1)=PP(:,i)+M*(a0*Q(:,i)+a2*V(:,i)+a3*AA(:,i))+C*(a1*Q(:,i)+a4*V(:,i)+a5*AA(:,i)); ik=L\PP(:,i+1); ik=U\ik; Q(:,i+1)=L'\ik; AA(:,i+1)=a0*(Q(:,i+1)-Q(:,i))-a2*V(:,i)-a3*AA(:,i); V(:,i+1)=V(:,i)+a6*AA(:,i)+a7*AA(:,i+1); end |
回复此楼» 猜你喜欢
售SCI一区T0P文章，我:8.O.55.1.O.5.4,科目全,可+急
已经有5人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有5人回复
-
时间戳又变了8-15
已经有11人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有5人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有3人回复
-
售SCI一区T0P文章，我:8.O55.1.O.54,科目全,可十急
已经有8人回复
-
欢迎发来filecode的Mz6后的代码验证其规律
已经有57人回复
-
有时候，自然基金真的不能太认真 （我的申报经验）
已经有9人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有3人回复
找到一些相关的精华帖子，希望有用哦~
如何用matlab求解矩阵系数的二阶微分方程
已经有13人回复
- matlab求解动力学 方程中的参数
已经有13人回复
- 请教matlab反应动力学参数估计遇到的问题，谢谢
已经有15人回复
- 【求助】催化反应动力学matlab计算各基元反应的速率常数时，该如何避免较小量被忽略？
已经有3人回复
点击这里搜索更多相关资源科研从小木虫开始，人人为我，我为人人
