版塊導(dǎo)航: 正在加載中...

應(yīng)《網(wǎng)絡(luò)安全法》要求，自2017年10月1日起，未進(jìn)行實(shí)名認(rèn)證將不得使用互聯(lián)網(wǎng)跟帖服務(wù)。為保障您的帳號(hào)能夠正常使用，請(qǐng)盡快對(duì)帳號(hào)進(jìn)行手機(jī)號(hào)驗(yàn)證，感謝您的理解與支持！

24小時(shí)熱門版塊排行榜

返回列表

caoxinchd

新蟲(chóng) (初入文壇)

應(yīng)助: 0 (幼兒園)
金幣: 17
帖子: 1
在線: 3.2小時(shí)
蟲(chóng)號(hào): 1342935
注冊(cè): 2011-07-11

[求助] 用Newmark方法計(jì)算系統(tǒng)的動(dòng)力學(xué)響應(yīng)的matlab程序

請(qǐng)大家?guī)兔纯催@個(gè)程序有什么問(wèn)題？用Newmark方法計(jì)算系統(tǒng)的動(dòng)力學(xué)響應(yīng)，結(jié)果大的驚人。
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);% 整體質(zhì)量矩陣
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系數(shù)
dt=0.01;
nt=20;%計(jì)算相應(yīng)時(shí)間
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;
%系數(shù)定義完
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

回復(fù)此樓

» 猜你喜歡

材料與化工（0856）304求 B區(qū) 調(diào)劑已經(jīng)有3人回復(fù)
268求調(diào)劑已經(jīng)有8人回復(fù)
265求調(diào)劑已經(jīng)有14人回復(fù)
350 本科985求調(diào)劑，求老登收留已經(jīng)有3人回復(fù)
【考研調(diào)劑】化學(xué)專業(yè) 281分，一志愿四川大學(xué)，誠(chéng)心求調(diào)劑已經(jīng)有9人回復(fù)
初試 317 已經(jīng)有3人回復(fù)
279分求調(diào)劑一志愿211 已經(jīng)有15人回復(fù)
286求調(diào)劑已經(jīng)有4人回復(fù)
279求調(diào)劑已經(jīng)有4人回復(fù)
求調(diào)劑已經(jīng)有3人回復(fù)

1樓 2013-05-29 11:54:44

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

智能機(jī)器人

Robot (super robot)

我們都愛(ài)小木蟲(chóng)

找到一些相關(guān)的精華帖子，希望有用哦~

如何用matlab求解矩陣系數(shù)的二階微分方程已經(jīng)有13人回復(fù)
matlab求解動(dòng)力學(xué) 方程中的參數(shù) 已經(jīng)有13人回復(fù)
請(qǐng)教matlab反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)遇到的問(wèn)題，謝謝已經(jīng)有15人回復(fù)
【求助】催化反應(yīng)動(dòng)力學(xué)matlab計(jì)算各基元反應(yīng)的速率常數(shù)時(shí)，該如何避免較小量被忽略？已經(jīng)有3人回復(fù)

點(diǎn)擊這里搜索更多相關(guān)資源

科研從小木蟲(chóng)開(kāi)始，人人為我，我為人人

相關(guān)版塊跳轉(zhuǎn) 我要訂閱樓主 caoxinchd 的主題更新