this is my code for one shaft 0 disk. but I don't know how to plot mode shape I can't find b.
nel1=14;
nnel=2; %
ndof=4; %
phase1=0; %
nb1=2; %
snode1=nel1+1; %
nnode1=(nel1+1); %
bpo1=[0,1.4]; %
edof=nnel*ndof; %
sdof=snode1*ndof; %
L1=1.4; %
Es1=2e11; % 
v1=0.3; % 
Ros1=7800; % 
Gs1=Es1/(2*(1+v1)); % 
G1=Gs1;
rr1=zeros(2,nel1); % 
rr1(1,:)=0.05; % 
rr1(2,:)=0; % 
MM1=zeros(sdof,sdof); % initialization of M system matrix
CC1=zeros(sdof,sdof); % initialization of C system matrix
KK1=zeros(sdof,sdof); % initialization of K system matrix
GG1=zeros(sdof,sdof)
gcoord1=zeros(nnode1,1);
for i=1:nnode1;
gcoord1(i)=gcoord1(i)+(i-1)/nel1*L1;
end
for i=1:nel1;
for j=1:nnel;
nodes1(i,j)=i+j-1;
end
end
for iel=1:nel1;
nl1=nodes1(iel,1); nr1=nodes1(iel,2);
L21=gcoord1(nr1)-gcoord1(nl1);
L11=gcoord1(nl1);
rso1=rr1(1,iel);
rsi1=rr1(2,iel);
m1=rsi1/rso1;
kpran1=(6*(1+v1)*(1+m1^2)^2)/((7+6*v1)*(1+m1^2)^2+(20+12*v1)*m1^2);
o1=(12*Es1*(1/4*pi*(rso1^4-rsi1^4)))/(kpran1*Gs1*(pi*(rso1^2-rsi1^2))*L21^2);
As1=pi*(rso1^2-rsi1^2);
Isd1=1/4*pi*(rso1^4-rsi1^4);
Isp1=1/2*pi*(rso1^4-rsi1^4);
[m1,k1,g1]=residualMCK1(o1,L1,Ros1,L11,L21,As1,Isp1,Isd1,Es1,kpran1,Gs1);
mm1=m1;
kk1=k1;
gg1=g1;
for i=1:edof
start1 = (iel-1)*(nnel-1)*ndof;
ii=start1+i;
for j=1:edof
jj=start1+j;
MM1(ii,jj)=MM1(ii,jj)+double(mm1(i,j));
KK1(ii,jj)=KK1(ii,jj)+double(kk1(i,j));
GG1(ii,jj)=GG1(ii,jj)+double(gg1(i,j));
end
end
end
KKB1=zeros(sdof,sdof);
KKB2=zeros(sdof,sdof);
CCB1=zeros(sdof,sdof);
CCB2=zeros(sdof,sdof);
Cyy1=0;
Cyy2=0;
Cyz1=0;
Czy1=0;
Czz1=0;
Czz2=0;
Kyy1=1e12;
Kyy2=1e12;
Kyz1=0;
Kzy1=0;
Kzz1=1e12;
Kzz2=1e12;
Kb1=[Kyy1 Kyz1;Kzy1 Kzz1];
Cb1=[Cyy1 Cyz1;Czy1 Czz1];
Kb2=[Kyy2 Kyz1;Kzy1 Kzz2];
Cb2=[Cyy2 Cyz1;Czy1 Czz2];
for ib1=1:nb1;
start1=bpo1(ib1)/L1*nel1*ndof;
for i=1:2;
ii=start1+i;
for j=1:2;
jj=start1+j;
KKB1(ii,jj)=KKB1(ii,jj)+Kb1(i,j); % combine bearing to global
CCB1(ii,jj)=CCB1(ii,jj)+Cb1(i,j);
KKB2(ii,jj)=KKB2(ii,jj)+Kb2(i,j); % combine bearing to global
CCB2(ii,jj)=CCB2(ii,jj)+Cb2(i,j);
end
end
end
MMMM=MM1;
CCCC1=CC1+CCB1;
CCCC2=CC1+CCB2;
KKKK1=KK1+KKB1;
KKKK2=KK1+KKB2;
GGGG=GG1;
o=zeros(sdof,sdof);
n=0;
XXYY=zeros(1,1);
Oma1=0;
n=n+1;
XXYY(n)=Oma1;%*60/2/pi;
E=[o,-MMMM;MMMM,(CCCC1)+(Oma1*GGGG)];
F=[MMMM,o;o,KKKK1];
EF=-inv(E)*F;
w=sort(abs(imag(eig(EF))));
wwwa=abs(imag(eig(EF)));
[C,D]=eig(EF);
K=imag(C);
R=real(C);
ww1(:,n)=w;%*60/2/pi;
A0=imag(D);
for k=1:sdof*2
A1(k,1)=A0(k,k);
end
for i=1:2:sdof*2-1
QQQ(:,i)=C(:,i)*D(i,i)+C(:,i+1)*D(i+1,i+1);
end
QQQ(:,sdof*2)=C(:,sdof*2)*D(sdof*2,sdof*2);
QQQV1=QQQ(sdof+1:ndof:sdof*2,:);
QQQW1=QQQ(sdof+2:ndof:sdof*2,:);
QQQB1=QQQ(sdof+3:ndof:sdof*2,:);
QQQR1=QQQ(sdof+4:ndof:sdof*2,:);
QQQV2=QQQ(sdof*1.5+1:ndof:sdof*2,:);
QQQW2=QQQ(sdof*1.5+2:ndof:sdof*2,:);
QQQU2=QQQ(sdof*1.5+3:ndof:sdof*2,:);
QQQB2=QQQ(sdof*1.5+4:ndof:sdof*2,:);
QQQR2=QQQ(sdof*1.5+5:ndof:sdof*2,:);
QQQA2=QQQ(sdof*1.5+6:ndof:sdof*2,:);
% plot mode shape % figure(1) %coupled b=117; plot(1:snode1,QQQV1(:,b),1:snode1,QQQW1(:,b),'linewidth',0.5); title('Mode 1 = 633.65 (rad/s)') legend('Driving-V','Driving-W') ;
figure(5)
%coupled
b=31;
plot(1:snode1,QQQV1(:,b),1:snode1,QQQW1(:,b),'linewidth',0.5);
title('Mode 2 = 2550 (rad/s)')
legend('Driving-V','Driving-W') ;

1 Kommentar

KSSV
KSSV am 12 Jul. 2017
We cannot run your code...few functions are missing.

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Kategorien

Mehr zu Develop Apps Using App Designer finden Sie in Hilfe-Center und File Exchange

Tags

Gefragt:

am 11 Jul. 2017

Kommentiert:

am 12 Jul. 2017

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by