Bardzo dziękuję za pomoc, bez metody newmarka raczej nie dałbym rady ;)
Poniżej cały kod, może komuś się przyda , rysowanie zrobione brute forcem - nie miałem już ochoty bawić się w pisanie pętli.
Pozdrawiam i życzę miłego wypoczynku przez resztę weekendu majowego ;)
l=4;
ro=7.85*10^3;
E=2.06*10^11;
%wprowadzanie wymuszen i war. początkowych
x=input('podaj wartość wychylenia środkowego bloczka od 0 do 1: ');
while x<0
disp(' !!! wychylenie musi być wartością dodatnią !!!');
while x>=1
disp(' !!! wychylenie musi być wartością mniejszą lub równą 1 !!!');
x=input('podaj wartość wychylenia środkowego bloczka od 0 do 1: ');
end
x=input('podaj wartość wychylenia środkowego bloczka od 0 do 1: ');
end
q=input('podaj wartość siły obciążającej środkowy bloczek (od -100 do 100): ');
while q>=100.1
disp(' !!! siła musi być wartością mniejszą lub równą 100 !!!');
q=input('podaj wartość siły obciążającej środkowy bloczek (od -100 do 100): ');
while q<=-100
disp(' !!! siła musi być wartością większą lub równą -100 !!!');
q=input('podaj wartość siły obciążającej środkowy bloczek (od -100 do 100): ');
end
end
%geometria el. skonczonego
L=l/5;
R=0.01;
R1=0.01;
A=3.14*R^2;
A1=3.14*R1^2;
L1=1.8;
J3=(3.14*R^4)/4;
J1=(3.14*R1^4)/4;
%boczne belki
MP=ro*L*[(13*A)/35+(6*J3)/(5*L^2) (11*A*L)/210+J3/(10*L) (9*A)/70-(6*J3)/(5*L^2) (-13*A*L)/420+J3/(10*L);
(11*A*L)/210+J3/(10*L) (A*L^2)/105+(2*J3)/15 (13*A*L)/420-J3/(10*L) (-A*L^2)/140-J3/30;
(9*A)/70-(6*J3)/(5*L^2) (13*A*L)/420-J3/(10*L) (13*A)/35+(6*J3)/(5*L^2) (-11*A*L)/210-J3/(10*L);
(-13*A*L)/420+J3/(10*L) (-A*L^2)/140-J3/30 (-11*A*L)/210-J3/(10*L) (A*L^2)/105+(2*J3)/15]
KP=(E/L)*[(12*J3)/(L^2) (6*J3)/L -(12*J3)/(L^2) (6*J3)/L;
(6*J3)/L 4*J3 -(6*J3)/L 2*J3;
-(12*J3)/(L^2) -(6*J3)/L (12*J3)/(L^2) -(6*J3)/L;
(6*J3)/L 2*J3 -(6*J3)/L 4*J3]
%środkowa belka
MP1=ro*L1*[(13*A1)/35+(6*J1)/(5*L1^2) (11*A1*L1)/210+J1/(10*L1) (9*A1)/70-(6*J1)/(5*L1^2) (-13*A1*L1)/420+J1/(10*L1);
(11*A1*L1)/210+J1/(10*L1) (A1*L1^2)/105+(2*J1)/15 (13*A1*L1)/420-J1/(10*L1) (-A1*L1^2)/140-J1/30;
(9*A1)/70-(6*J1)/(5*L1^2) (13*A1*L1)/420-J1/(10*L1) (13*A1)/35+(6*J1)/(5*L1^2) (-11*A1*L1)/210-J1/(10*L1);
(-13*A1*L1)/420+J1/(10*L1) (-A1*L1^2)/140-J1/30 (-11*A1*L1)/210-J1/(10*L1) (A1*L1^2)/105+(2*J1)/15]
KP1=(E/L1)*[(12*J1)/(L1^2) (6*J1)/L1 -(12*J1)/(L1^2) (6*J1)/L1;
(6*J1)/L1 4*J1 -(6*J1)/L1 2*J1;
-(12*J1)/(L1^2) -(6*J1)/L1 (12*J1)/(L1^2) -(6*J1)/L1;
(6*J1)/L1 2*J1 -(6*J1)/L1 4*J1]
M=zeros(24,24);
K=zeros(24,24);
p=0;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=2;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=4;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=6;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=8;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=10;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP1(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP1(iw,ik)+K(iw+p,ik+p);
end
end
p=12;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=14;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=16;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=18;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
p=20;
for iw=1:4
for ik=1:4
M(iw+p,ik+p)=MP(iw,ik)+M(iw+p,ik+p);
K(iw+p,ik+p)=KP(iw,ik)+K(iw+p,ik+p);
end
end
K(:,[1 2 23 24])=0;
K([1 2 23 24],:)=0;
K(1,1)=1;
K(2,2)=1;
K(23,23)=1;
K(24,24)=1;
M(:,[1 2 23 24])=0;
M([1 2 23 24],:)=0;
M(1,1)=1;
M(2,2)=1;
M(23,23)=1;
M(24,24)=1;
Q=zeros(1,24)';
Q(12,1)=q/2;
Q(14,1)=q/2;
U=inv(K)*Q;
omega2 = eig(K,M);
fmes = sqrt(omega2(2))/(2*pi);
omega = ( pi/(2*(L+1/2*L1)) ) * sqrt(E/ro);
f = omega/(2*pi) ;
% -----------------------------------------------------------
[r,k]=size(K);
t=5;
% time step
dt=0.1;
% Przemieszczenie
u0=zeros(r,1);
u0(11,1)=x;
u0(13,1)=x;
% Prędkość
udot0=zeros(r,1);
% Przyspieszenie
u2dot0=zeros(r,1);
% Wymuszenie (Siła)
R0=zeros(r,0);
R0(11,1)=q/2;
R0(13,1)=q/2;
% parametry
beta=0.5;
alfa=0.25*(0.5+beta)^2;
% Stałe
a0=1/(alfa*dt^2);
a1=beta/(alfa*dt);
a2=1/(alfa*dt);
a3=1/(2*alfa)-1;
a4=beta/alfa-1;
a5=dt/2*(beta/alfa-2);
a6=dt*(1-beta);
a7=beta*dt;
% Efektywna macierz sztywności
K_e=K+a0*M;
u=horzcat(u0,zeros(r,length(dt:dt:t)-1));
udot=horzcat(udot0,zeros(r,length(dt:dt:t)-1));
u2dot=horzcat(u2dot0,zeros(r,length(dt:dt:t)-1));
R=horzcat(R0,zeros(r,length(dt:dt:t)-1));
for i=1:1:length(dt:dt:t)-1;
% effective load
R_e=R(:,i)+M*(a0*u(:,i)+a2*udot(:,i)+a3*u2dot(:,i));
u(:,i+1)=inv(K_e)*R_e;
u2dot(:,i+1)=a0*(u(:,i+1)-u(:,i))-a2*udot(:,i)-a3*u2dot(:,i);
udot(:,i+1)=udot(:,i)+a6*u2dot(:,i)+a7*u2dot(:,i+1);
end
clear a0 a1 a2 a3 a4 a5 a6 a7
clear K_e u0 udot0 u2dot0 R0
clear alfa beta
u2=u(3,:)
rot2=u(4,:)
u3=u(5,:)
rot3=u(6,:)
u4=u(7,:)
rot4=u(8,:)
u5=u(9,:)
rot5=u(10,:)
u6=u(11,:)
rot6=u(12,:)
u7=u(13,:)
rot7=u(14,:)
u8=u(15,:)
rot8=u(16,:)
u9=u(17,:)
rot9=u(18,:)
u10=u(19,:)
rot10=u(20,:)
u11=u(21,:)
rot11=u(22,:)
jj_end=max(size(u));
for jj=1:1:jj_end
hold off
plot(0,0)
hold on
%el1
plot([0.2,1-(tan(rot2(jj))*0.5)],[2.8,2.8+u2(jj)])
plot([0.2,1+(tan(rot2(jj))*0.5)],[3.2,3.2+u2(jj)])
plot([1-(tan(rot2(jj))*0.5),1+(tan(rot2(jj))*0.5)],[2.8+u2(jj),3.2+u2(jj)])
%el2
plot([1-tan(rot2(jj))*0.5,1.8-(tan(rot3(jj))*0.5)],[2.8+u2(jj),2.8+u3(jj)])
plot([1+tan(rot2(jj))*0.5,1.8+(tan(rot3(jj))*0.5)],[3.2+u2(jj),3.2+u3(jj)])
plot([1.8-tan(rot3(jj))*0.5,1.8+tan(rot3(jj))*0.5],[2.8+u3(jj),3.2+u3(jj)])
%el3
plot([1.8-tan(rot3(jj)*0.5),2.6-(tan(rot4(jj))*0.5)],[2.8+u3(jj),2.8+u4(jj)])
plot([1.8+tan(rot3(jj)*0.5),2.6+(tan(rot4(jj))*0.5)],[3.2+u3(jj),3.2+u4(jj)])
plot([2.6-tan(rot4(jj)*0.5)*0.5,2.6+tan(rot4(jj))*0.5],[2.8+u4(jj),3.2+u4(jj)])
%el4
plot([2.6-tan(rot4(jj))*0.5,3.4-(tan(rot5(jj))*0.5)],[2.8+u4(jj),2.8+u5(jj)])
plot([2.6+tan(rot4(jj))*0.5,3.4+(tan(rot5(jj))*0.5)],[3.2+u4(jj),3.2+u5(jj)])
plot([3.4-tan(rot5(jj))*0.5,3.4+tan(rot5(jj))*0.5],[2.8+u5(jj),3.2+u5(jj)])
%el5
plot([3.4-tan(rot5(jj))*0.5,4.2],[2.8+u5(jj),2.8+u6(jj)])
plot([3.4+tan(rot5(jj))*0.5,4.2],[3.2+u5(jj),3.2+u6(jj)])
%el6 (uproszczony) brak rotacji ścian el 5 i 7 przyczepione na
%stałe do ścainek bocznych.
plot([4.2,6.2,6.2,4.2,4.2],[2+u6(jj),2+u6(jj),4+u7(jj),4+u7(jj),2+u6(jj)])
%el7
plot([7-tan(rot8(jj))*0.5,6.2],[2.8+u8(jj),2.8+u7(jj)])
plot([7+tan(rot8(jj))*0.5,6.2],[3.2+u8(jj),3.2+u7(jj)])
plot([7-tan(rot8(jj))*0.5,7+tan(rot8(jj))*0.5],[2.8+u8(jj),3.2+u8(jj)])
%el8
plot([7-tan(rot8(jj))*0.5,7.8-(tan(rot9(jj))*0.5)],[2.8+u8(jj),2.8+u9(jj)])
plot([7+tan(rot8(jj))*0.5,7.8+(tan(rot9(jj))*0.5)],[3.2+u8(jj),3.2+u9(jj)])
plot([7.8-tan(rot9(jj))*0.5,7.8+tan(rot9(jj))*0.5],[2.8+u9(jj),3.2+u9(jj)])
%el9
plot([7.8-tan(rot9(jj))*0.5,8.6-(tan(rot10(jj))*0.5)],[2.8+u9(jj),2.8+u10(jj)])
plot([7.8+tan(rot9(jj))*0.5,8.6+(tan(rot10(jj))*0.5)],[3.2+u9(jj),3.2+u10(jj)])
plot([8.6-tan(rot10(jj))*0.5,8.6+tan(rot10(jj))*0.5],[2.8+u10(jj),3.2+u10(jj)])
%el10
plot([8.6-tan(rot10(jj))*0.5,9.4-(tan(rot11(jj))*0.5)],[2.8+u10(jj),2.8+u11(jj)])
plot([8.6+tan(rot10(jj))*0.5,9.4+(tan(rot11(jj))*0.5)],[3.2+u10(jj),3.2+u11(jj)])
plot([9.4-tan(rot10(jj))*0.5,9.4+tan(rot10(jj))*0.5],[2.8+u11(jj),3.2+u11(jj)])
%el11
plot([9.4-tan(rot11(jj))*0.5,10.2],[2.8+u11(jj),2.8])
plot([9.4+tan(rot11(jj))*0.5,10.2],[3.2+u11(jj),3.2])
%sciana
plot(0,0,5,5, [0.2,0.2],[2.3,3.6])
plot(0,0,5,5, [10.2,10.2],[2.3,3.6])
%kreskowanie
for ii=9:12
plot([0.2 0.1], [4/15*ii 4/15*ii+0.1])
end
for ii=9:12
plot([10.3 10.2], [4/15*ii+0.1 4/15*ii])
end
pause(0.05)
end
disp('Siła wymuszenia=');
q
disp('Przesunięcie środkowego bloczka=');
x
disp('Częstość własna układu');
omega
|
Użytkownicy Online
ew thread" po prawej stronie.
