System Differentialgleichunge < Matlab < Mathe-Software < Mathe < Vorhilfe
|
Hallo
ich hab eine Bewegungsgleichung mit 6 Freiheitsgraden, dementsprechend natürlich ein Gleichungssystem aus 6 Differentialgleichungen.
form: Mq´´+Dq´+Kq = f(t)
für MAtlab in der Form : q´´= Mi(f(t)-Dq´-kq) Mi = invertierte Massenmatrix
Ich habe MatLab soweit gebracht, dass er mir die Gleichungen löst, allerdings steckt irgendwo noch (mindestens) ein Fehler.
das hier ist der Quellcode zur Lösung des Systems:
y0 = [0 ; 0 ; 0.001 ;0 ; 0 ; 0 ; 0 ;0 ; 0; 0; 0 ;0] % Anfangswerte
tint = [0 1] % Integrationsintervall
[T,Y] = ode45 (@Bewegungsgleichung, tint, y0)
sollte soweit stimmen.
meine Funktion heißt Bewegungsgleichung. im moment ist das system ungedämpft, es gibt keine zusätzlichen kräfte (ensprechende Teile im Quellcode mit % abgetrennt) also in der Form: q´´= Mi*K*q
mein Problem ist, dass die Anfangswerte nicht angenommen werden und yp(3) kein spalten, sondern zeilenvektor wird. Wäre super, wenn mir jemand helfen könnte.
das ist die Funktion:
function yp = Bewegungsgleichung (t,y)
yp = zeros(6,1) %macht yp zu column Vector
% Variablen
E1 = 200000 % E-Modul Stange
E2 = 200000 % E-Modul Welle
d1 = 0.01 % Durchmesser der Welle
d2 = 0.02 % Durchmesser der Stange
l = 0.5 % effektive Länge Stange, Welle
l1 = 0.135 % Abstand Kraftangriffspunk/Einspannung
l2 = l-l1
m1 = 23 % Masse Oberplatte + Motor
m2 = [mm] pi*(d1/2)^2*l2*7850 [/mm] % Massse an q2
m3 = 7.52 % Masse Shaker + Halter
m4 = [mm] pi*(d2/2)^2*l*7850 [/mm] % MAsse an q4
m6 = [mm] pi*(d1/2)^2*l1*7850 [/mm] % Masse an q6
I1 = (pi * [mm] (d2/2)^2)/4 [/mm] % Flächenträgheitsmoment Stange
I2 = (pi * [mm] (d1/2)^2)/4 [/mm] % Flächenträgheitsmoment Welle
% Komponenten der Steifigkeitsmatrix
k11 = [mm] 24/l^3*(2*E1*I1+E2*I2) [/mm] % aus q1 = [1 0 0 0 0 0]
k12 = [mm] -24*E2*I2/l^3
[/mm]
k13 = 0
k14 = [mm] -24*E1*I1/l^3
[/mm]
k15 = [mm] -24*E1*I1/l^3
[/mm]
k16 = 0
k21 = k12 % aus q2 = [0 1 0 0 0 0]
k22 = [mm] 768*E2*I2/(7/l2^3)
[/mm]
k23 = 0
k24 = 0
k25 = 0
k26 = [mm] k22*(l2/l)^2*(1+l1/(2*l))
[/mm]
k31 = k13
k32 = k23
k33 = [mm] 4*E1*I1*l^2/(l1^2*l2^3*(1+l1/(3*l)))
[/mm]
k34 = [mm] k33*(l2/l)^2*(1+l1/(2*l))
[/mm]
k35 = 0
k36 = 0
k41 = k14
k42 = k24
k43 = k34
k44 = [mm] 768*E2*I2/(7*l2^3)
[/mm]
k45 = 0
k46 = 0
k51 = k15
k52 = k25
k53 = k35
k54 = k45
k55 = [mm] 192*E1*I1/l^3
[/mm]
k56 = 0
k61 = k16
k62 = k26
k63 = k36
k64 = k46
k65 = k56
k66 = [mm] 4*E2*I2*l^2/(l1^2*l2^3*(1+l1/(3*l)))
[/mm]
%Steifigkeitsmatrix
K = [k11 k21 k31 k41 k51 k61; k12 k22 k32 k42 k52 k62; k13 k23 k33 k43 k53 k63; k14 k24 k34 k44 k54 k64
; k15 k25 k35 k45 k55 k65 ; k16 k26 k36 k46 k56 k66]
%Elemete der Dämpfungsmatrix
d11 = 1 % aus q1 = [1 0 0 0 0 0]
d12 = 1
d13 = 0
d14 = 1
d15 = 1
d16 = 0
d21 = d12 % aus q2 = [0 1 0 0 0 0]
d22 = 1
d23 = 0
d24 = 0
d25 = 0
d26 = 1
d31 = d13
d32 = d23
d33 = 1
d34 = 1
d35 = 0
d36 = 0
d41 = d14
d42 = d24
d43 = d34
d44 = 1
d45 = 0
d46 = 0
d51 = d15
d52 = d25
d53 = d35
d54 = d45
d55 = 1
d56 = 0
d61 = d16
d62 = d26
d63 = d36
d64 = d46
d65 = d56
d66 = 1
%Dämpfungsmatrix
D = [d11 d21 d31 d41 d51 d61; d12 d22 d32 d42 d52 d62; d13 d23 d33 d43 d53 d63; d14 d24 d34 d44 d54 d64
; d15 d25 d35 d45 d55 d65 ; d16 d26 d36 d46 d56 d66]
%Massenmatrix
M = [m1 0 0 0 0 0; 0 m2 0 0 0 0; 0 0 m3 0 0 0; 0 0 0 m4 0 0; 0 0 0 0 m4 0; 0 0 0 0 0 m6]
Mi = inv(M)
% Matrix F = Mi*K
F = Mi*K
FT = transpose (F)
% Matrix G = Mi*D
G = Mi*D
GT = transpose (G)
% Shakerkraft
% x = pi:0.01:pi
% f = 9.1*sin(x)
%Gleichungssystem übeführt in System erster ordnung
yp(1) = y(7)
yp(2) = y(8)
YP(3) = y(9)
yp(4) = y(10)
yp(5) = y(11)
yp(6) = y(12)
yp(7) = -((transpose (FT*[1;0;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (FT*[1;0;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(8) = -((transpose (FT*[0;1;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (FT*[0;1;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(9) = -((transpose (FT*[0;0;1;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (FT*[0;0;1;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %+f %+gottverdammte scheißkraft
yp(10) = -((transpose (FT*[0;0;0;1;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (FT*[0;0;0;1;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(11) = -((transpose (FT*[0;0;0;0;1;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (FT*[0;0;0;0;1;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(12) = -((transpose (FT*[0;0;0;0;0;1]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (FT*[0;0;0;0;0;1]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %-f %-gottverdammte scheißkraft
end
ch habe diese Frage auch in folgenden Foren auf anderen Internetseiten gestellt:
http://www.maschinenbauer-forum.de/thread.php?threadid=9202&sid=
|
|
| |
|
Hi,
ich krieg das schon irgendwie gelöst, und die Lösung hängt auch irgendwie
von den Anfangsbedingungen ab.
Aber Du solltest überlegen, ob Du mit ein paar ; nicht hier und da die Ausgabe
unterdrücken möchtest, ob Du bei jedem Funktionsaufruf von Bewegungs-
gleichungen.m immer transponieren musst. Evtl ginge das mit globalen
Variablen schneller.
Ausserden ist eines von den yp's gross geschrieben.
Und ich glaube in der Hilfe zu ode* steht, dass man die Startbedingungen
als yp0=[ x x x x] eingibt, also als Zeilenvektor. Das ist aber glaube ich egal.
mfg
nschlange
|
|
|
| |
|
ich habe nochmal gekuckt, die anfangsbedingungen werden so eingegeben, wie ich es habe.
und du hast recht, irgendwie verändrt sich das ganze schon mit den startbedingungen, allerdings nicht so, wie es sein sollte. zum zeitpunkt t=0 sollte ja der jeweilige startwert auch angezeigt und berücksichtigt werden. passiert aber leider nicht. insbesondere y(3) bleibt meist konstant auf null, auch wenn ich einen Anfangswert eingebe.
zu den ; ,sorry die habe ich in der neusten version gesetzt, nur in der gepostetten hatte ich es noch nicht, einfach um alle ergebnisse mal auf die plausibilität prüfen zu können. hier nochmal mit ; :
%Bewegungsgleichung
function yp = Bewegungsgleichung (t,y)
yp = zeros(6,1) %macht yp zu column Vector
% Variablen
E1 = 200000 % E-Modul Stange
E2 = 200000 % E-Modul Welle
d1 = 0.01 % Durchmesser der Welle
d2 = 0.02 % Durchmesser der Stange
l = 0.5 % effektive Länge Stange, Welle
l1 = 0.135 % Abstand Kraftangriffspunk/Einspannung
l2 = l-l1
m1 = 23 % Masse Oberplatte + Motor
m2 = [mm] pi*(d1/2)^2*l2*7850 [/mm] % Massse an q2
m3 = 7.52 % Masse Shaker + Halter
m4 = [mm] pi*(d2/2)^2*l*7850 [/mm] % MAsse an q4
m6 = [mm] pi*(d1/2)^2*l1*7850 [/mm] % Masse an q6
I1 = (pi * [mm] (d2/2)^2)/4 [/mm] % Flächenträgheitsmoment Stange
I2 = (pi * [mm] (d1/2)^2)/4 [/mm] % Flächenträgheitsmoment Welle
% Komponenten der Steifigkeitsmatrix
k11 = [mm] 24/l^3*(2*E1*I1+E2*I2) [/mm] ; % aus q1 = [1 0 0 0 0 0]
k12 = [mm] -24*E2*I2/l^3;
[/mm]
k13 = 0;
k14 = [mm] -24*E1*I1/l^3;
[/mm]
k15 = [mm] -24*E1*I1/l^3;
[/mm]
k16 = 0;
k21 = k12; % aus q2 = [0 1 0 0 0 0]
k22 = [mm] 768*E2*I2/(7/l2^3);
[/mm]
k23 = 0;
k24 = 0;
k25 = 0;
k26 = [mm] k22*(l2/l)^2*(1+l1/(2*l));
[/mm]
k31 = k13;
k32 = k23;
k33 = [mm] 4*E1*I1*l^2/(l1^2*l2^3*(1+l1/(3*l)));
[/mm]
k34 = [mm] k33*(l2/l)^2*(1+l1/(2*l));
[/mm]
k35 = 0;
k36 = 0;
k41 = k14;
k42 = k24;
k43 = k34;
k44 = [mm] 768*E2*I2/(7*l2^3);
[/mm]
k45 = 0;
k46 = 0;
k51 = k15;
k52 = k25;
k53 = k35;
k54 = k45;
k55 = [mm] 192*E1*I1/l^3;
[/mm]
k56 = 0;
k61 = k16;
k62 = k26;
k63 = k36;
k64 = k46;
k65 = k56;
k66 = [mm] 4*E2*I2*l^2/(l1^2*l2^3*(1+l1/(3*l)));
[/mm]
%Steifigkeitsmatrix
K = [k11 k21 k31 k41 k51 k61; k12 k22 k32 k42 k52 k62; k13 k23 k33 k43 k53 k63; k14 k24 k34 k44 k54 k64
; k15 k25 k35 k45 k55 k65 ; k16 k26 k36 k46 k56 k66]
%Elemete der Dämpfungsmatrix
d11 = 1; % aus q1 = [1 0 0 0 0 0]
d12 = 1;
d13 = 0;
d14 = 1;
d15 = 1;
d16 = 0;
d21 = d12; % aus q2 = [0 1 0 0 0 0]
d22 = 1;
d23 = 0;
d24 = 0;
d25 = 0;
d26 = 1;
d31 = d13;
d32 = d23;
d33 = 1;
d34 = 1;
d35 = 0;
d36 = 0;
d41 = d14;
d42 = d24;
d43 = d34;
d44 = 1;
d45 = 0;
d46 = 0;
d51 = d15;
d52 = d25;
d53 = d35;
d54 = d45;
d55 = 1;
d56 = 0;
d61 = d16;
d62 = d26;
d63 = d36;
d64 = d46;
d65 = d56;
d66 = 1;
%Dämpfungsmatrix
D = [d11 d21 d31 d41 d51 d61; d12 d22 d32 d42 d52 d62; d13 d23 d33 d43 d53 d63; d14 d24 d34 d44 d54 d64
; d15 d25 d35 d45 d55 d65 ; d16 d26 d36 d46 d56 d66]
%Massenmatrix
M = [m1 0 0 0 0 0; 0 m2 0 0 0 0; 0 0 m3 0 0 0; 0 0 0 m4 0 0; 0 0 0 0 m4 0; 0 0 0 0 0 m6]
Mi = inv(M) %invertierte Massenmatrix
% Matrix F = Mi*K
F = Mi*K;
FT = transpose (F);
% Matrix G = Mi*D
G = Mi*D;
GT = transpose (G);
% Shakerkraft
% x = -pi:0.01:pi
% f = 9.1*sin(x)
%Gleichungssystem übeführt in System erster ordnung
yp(1) = y(7)
yp(2) = y(8)
YP(3) = y(9)
yp(4) = y(10)
yp(5) = y(11)
yp(6) = y(12)
yp(7) = -((transpose (FT*[1;0;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (GT*[1;0;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(8) = -((transpose (FT*[0;1;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (GT*[0;1;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(9) = -((transpose (FT*[0;0;1;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (GT*[0;0;1;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %+f %+shakerkraft
yp(10) = -((transpose (FT*[0;0;0;1;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (GT*[0;0;0;1;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(11) = -((transpose (FT*[0;0;0;0;1;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (GT*[0;0;0;0;1;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(12) = -((transpose (FT*[0;0;0;0;0;1]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) % -((transpose (GT*[0;0;0;0;0;1]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %-f %-shakerkraft
end
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 12:12 Do 09.08.2007 | | Autor: | nschlange |
Ich weiss es jetzt nicht genau, aber es könnte daran liegen, dass du yp(3) als YP(3) geschrieben hast. Zu mindest in der copy-paste Version.
Probier doch mal das zu ändern.
|
|
|
| |
|
hey danke, das hat schonmal geholfen das problem mit dem zeilen/spaltenvektor bei yp(3) in den griff zu kriegen.
die anfangswerte stimmen aber leider immer noch nicht.
|
|
|
| |
|
Hallo
also, ich habe das große y korrigiert (danke für den tip, solche fehler sieht man einfach irgendwann nicht mehr), und jetzt das problem mit dem spaltem/zeilenvektor gelöst. Leider nimmt matlab immer noch nichtt die anfangswerte an. wer hat noch eine idee woran das liegebn könnte?
|
|
|
| |
|
Hi,
wähl mal als Integrationsintervall 0..0.2 oder so.
Die Werte werden so gross, Größenordnung [mm] 10^6, [/mm] so dass,
wenn Du alle Komponenten von Y in einem Diagramm darstellst,
keine Unterschiede ausmachen kannst.
Vielleicht ist ja auch die DGL falsch.
mfg
nschlange
|
|
|
| |
|
Hallo
also mein system funktioniert jetzt soweit. das heißt, es kommen für startwerte vernünftige lösungen raus.
jetzt würde ich aber gerne das sytem durch eine äußere kraft anregen, z.B durch einen sinus. mein problem ist, das meine variable im sinus sich mit den zeitschritten der integration verändert, ich also z. B von 0-2pi gehen will und matlab mir das in genau soviele zeitschritte einteilt, wie er für die lösung der dgl benutzt. kann ich die zeitschritte evtl festlegen? oder kann ich die länge des zeitvektors auslesen?
hier nochmal das programm:
y0 = [0.001 ; 0 ; 0.00 ; 0 ; 0 ; -0.001 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0] % Anfangswerte
tint = [0 200] % Integrationsintervall
[T,Y] = ode45 (@Bewegungsgleichung, tint, y0) %Integrationsbefehl
und die funktion Bewegungsgleichung: diese funktion funktioniert, außer, wie schon oben erwähnt die erregung durch die kraft
function yp = Bewegungsgleichung (t,y,f)
yp = zeros(6,1) %macht yp zu column Vector
% Variablen
E1 = 0.200000 ; % E-Modul Stange (N/m²)
E2 = 0.200000 ; % E-Modul Welle
d1 = 0.01 ; % Durchmesser der Welle
d2 = 0.02 ; % Durchmesser der Stange
l = 0.5 ; % effektive Länge Stange, Welle
l1 = 0.135 ; % Abstand Kraftangriffspunk/Einspannung
l2 = l-l1;
m1 = 23 ; % Masse Oberplatte + Motor
m2 = [mm] pi*(d1/2)^2*l2*7850; [/mm] % Massse an q2
m3 = 7.52 ; % Masse Shaker + Halter
m4 = [mm] pi*(d2/2)^2*l*7850 [/mm] ; % Masse an q4
m6 = [mm] pi*(d1/2)^2*l1*7850 [/mm] ; % Masse an q6
I1 = (pi * [mm] (d2/2)^2)/4 [/mm] ; % Flächenträgheitsmoment Stange
I2 = (pi * [mm] (d1/2)^2)/4 [/mm] ; % Flächenträgheitsmoment Welle
% Komponenten der Steifigkeitsmatrix
k11 = [mm] 96/l^3*(4*E1*I1+E2*I2) [/mm] ; % aus q1 = [1 0 0 0 0 0]
k12 = [mm] -96*E2*I2/l^3;
[/mm]
k13 = 0;
k14 = [mm] -96*E1*I1/l^3;
[/mm]
k15 = [mm] -96*E1*I1/l^3;
[/mm]
k16 = 0;
k21 = k12; % aus q2 = [0 1 0 0 0 0]
k22 = [mm] 3*E2*I2*(l2)^3/((l/2)^3*(l2-l/2)^3);
[/mm]
k23 = 0;
k24 = [mm] k22*((l2-l/2)/l2)^2*(1+l/l2);
[/mm]
k25 = k24;
k26 = -k24;
k31 = k13;
k32 = k23;
k33 = [mm] 3*E1*I1*(l/2)^3/((l1)^3*(l/2-l1)^3);
[/mm]
k34 = [mm] k33*(l2/l)^2*(1+l1/(2*l));
[/mm]
k35 = 0;
k36 = 0;
k41 = k14;
k42 = k24;
k43 = k34;
k44 = [mm] 3*E1*I1*l2^3/((l2-l/2)^3*(l/2)^3);
[/mm]
k45 = k42;
k46 = 0;
k51 = k15;
k52 = k25;
k53 = k35;
k54 = k45;
k55 = [mm] 192*E1*I1/l^3;
[/mm]
k56 = 0;
k61 = k16;
k62 = k26;
k63 = k36;
k64 = k46;
k65 = k56;
k66 = [mm] 3*E2*I2*(l/2)^3/((l/2-l1)^3*l1^3);
[/mm]
%Steifigkeitsmatrix
K = [k11 k21 k31 k41 k51 k61; k12 k22 k32 k42 k52 k62; k13 k23 k33 k43 k53 k63; k14 k24 k34 k44 k54 k64
; k15 k25 k35 k45 k55 k65 ; k16 k26 k36 k46 k56 k66]
%Elemete der Dämpfungsmatrix
d11 = 1; % aus q1 = [1 0 0 0 0 0]
d12 = 1;
d13 = 0;
d14 = 1;
d15 = 1;
d16 = 0;
d21 = d12; % aus q2 = [0 1 0 0 0 0]
d22 = 1;
d23 = 0;
d24 = 1;
d25 = 1;
d26 = 1;
d31 = d13;
d32 = d23;
d33 = 1;
d34 = 1;
d35 = 0;
d36 = 0;
d41 = d14;
d42 = d24;
d43 = d34;
d44 = 1;
d45 = 1;
d46 = 0;
d51 = d15;
d52 = d25;
d53 = d35;
d54 = d45;
d55 = 1;
d56 = 0;
d61 = d16;
d62 = d26;
d63 = d36;
d64 = d46;
d65 = d56;
d66 = 1;
%Dämpfungsmatrix
D = [d11 d21 d31 d41 d51 d61; d12 d22 d32 d42 d52 d62; d13 d23 d33 d43 d53 d63; d14 d24 d34 d44 d54 d64
; d15 d25 d35 d45 d55 d65 ; d16 d26 d36 d46 d56 d66]
%Massenmatrix
M = [m1 0 0 0 0 0; 0 m2 0 0 0 0; 0 0 m3 0 0 0; 0 0 0 m4 0 0; 0 0 0 0 m4 0; 0 0 0 0 0 m6]
Mi = inv(M) %invertierte Massenmatrix
% Matrix F = Mi*K
F = Mi*K;
FT = transpose (F);
% Matrix G = Mi*D
G = Mi*D;
GT = transpose (G);
% Shakerkraft
x = -pi:(2*pi)/length (T):pi
f = 9.1*sin (t/t*x)
%Gleichungssystem übeführt in System erster ordnung
yp(1) = y(7)
yp(2) = y(8)
yp(3) = y(9)
yp(4) = y(10)
yp(5) = y(11)
yp(6) = y(12)
yp(7) = -((transpose (FT*[1;0;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[1;0;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(8) = -((transpose (FT*[0;1;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[0;1;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(9) = -((transpose (FT*[0;0;1;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)])+f %-((transpose (GT*[0;0;1;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %+f %+shakerkraft
yp(10) = -((transpose (FT*[0;0;0;1;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[0;0;0;1;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(11) = -((transpose (FT*[0;0;0;0;1;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[0;0;0;0;1;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(12) = -((transpose (FT*[0;0;0;0;0;1]))*[y(1);y(2);y(3);y(4);y(5);y(6)])-f %-((transpose (GT*[0;0;0;0;0;1]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %-f %-shakerkraft
end
|
|
|
| |
|
Hi Hendrik!!
Also dein Zeitvektor steckt doch in deinem "T" drin!! Das ist der Vektor, der von deiner DGL angelegt wird!
--------------------------------------------------------------
Ansonsten:
Schau dir mal deine Shakerkraft an:
% Shakerkraft
x = -pi:(2*pi)/length (T):pi
f = 9.1*sin (t/t*x)
--------------------------------------------------------------
Meiner Meinung nach ist f nur abhängig von x, soll heißen: (t/t*x)=x
Klammer falsch gesetzt?????
Viel Spaß noch! 
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 15:07 Mi 15.08.2007 | | Autor: | NewtonsLaw |
Als fehlerhaft gekennzeichnet, weil ich mir nicht sicher bin!
|
|
|
| |
|
hallo
ja, den vekot T kann ich bekommen, auch seine länge. Das problem ist, das dieser Vektor erst existiert, wenn die Integration gelaufen ist, ich also nicht wirklich schon vorher damitr arbeiten kann.
|
|
|
| |
|
Ich formuliere die Frage einfach nochmal anders. Ich habe die Lösungen der Homogenen Gleichung, (also mq´´+dq´+kq = 0).
mich interessiert aber die Lösung der inhomogenen Gleichung, also mit einer Kraft mit verschiedenen Frequenzen auf der rechten Seite
(also mq´´+dq´+kq = f(t))
wer hat eine idee, wie das funktionieren könnte?
Hallo
hab die lösung gefunden.
hier für interessierte einfach nochmal das gesamte programm.
mfg
Hendrik
Der Integrationsaufruf
% Berechnung des Gesamtsystems Demonstrator
%Lösung der Diffenretialgleichung Mq´´+ Dq´+ kq = f(t)
y0 = [0.000 ; 0 ; 0.00 ; 0 ; 0 ; -0.00 ; 0 ; 0 ; 0 ; 0 ; 0 ; 0] % Anfangswerte
tint = linspace(0, 100, 200) % [0 , 100]Integrationsintervall
[T,Y] = ode45 (@Bewegungsgleichung, tint, y0) %Integrationsbefehl
%Visualisierung
plot(T,Y(:,1),'-',T,Y(:,2),'-',T,Y(:,3),'-.',T,Y(:,4),'-',T,Y(:,5),'-',T,Y(:,6),'--')
title('Ausschlag über der Zeit');
ylabel('Ausschlag X');
xlabel('Zeit T');
[mm] legend('x_1','x_2','x_3','x_4','x_5','x_6')
[/mm]
und die funktion
function yp = Bewegungsgleichung (t,y)
yp = zeros(6,1) %macht yp zu column Vector
% Variablen
E1 = 0.200000 ; % E-Modul Stange (N/m²)
E2 = 0.200000 ; % E-Modul Welle
d1 = 0.01 ; % Durchmesser der Welle
d2 = 0.02 ; % Durchmesser der Stange
l = 0.5 ; % effektive Länge Stange, Welle
l1 = 0.135 ; % Abstand Kraftangriffspunk/Einspannung
l2 = l-l1;
m1 = 23 ; % Masse Oberplatte + Motor
m2 = [mm] pi*(d1/2)^2*l2*7850; [/mm] % Massse an q2
m3 = 7.52 ; % Masse Shaker + Halter
m4 = [mm] pi*(d2/2)^2*l*7850 [/mm] ; % Masse an q4
m6 = [mm] pi*(d1/2)^2*l1*7850 [/mm] ; % Masse an q6
I1 = (pi * [mm] (d2/2)^2)/4 [/mm] ; % Flächenträgheitsmoment Stange
I2 = (pi * [mm] (d1/2)^2)/4 [/mm] ; % Flächenträgheitsmoment Welle
% Komponenten der Steifigkeitsmatrix
k11 = [mm] 96/l^3*(4*E1*I1+E2*I2) [/mm] ; % aus q1 = [1 0 0 0 0 0]
k12 = [mm] -96*E2*I2/l^3;
[/mm]
k13 = 0;
k14 = [mm] -96*E1*I1/l^3;
[/mm]
k15 = [mm] -96*E1*I1/l^3;
[/mm]
k16 = 0;
k21 = k12; % aus q2 = [0 1 0 0 0 0]
k22 = [mm] 3*E2*I2*(l2)^3/((l/2)^3*(l2-l/2)^3);
[/mm]
k23 = 0;
k24 = [mm] k22*((l2-l/2)/l2)^2*(1+l/l2);
[/mm]
k25 = k24;
k26 = -k24;
k31 = k13;
k32 = k23;
k33 = [mm] 3*E1*I1*(l/2)^3/((l1)^3*(l/2-l1)^3);
[/mm]
k34 = [mm] k33*(l2/l)^2*(1+l1/(2*l));
[/mm]
k35 = 0;
k36 = 0;
k41 = k14;
k42 = k24;
k43 = k34;
k44 = [mm] 3*E1*I1*l2^3/((l2-l/2)^3*(l/2)^3);
[/mm]
k45 = k42;
k46 = 0;
k51 = k15;
k52 = k25;
k53 = k35;
k54 = k45;
k55 = [mm] 192*E1*I1/l^3;
[/mm]
k56 = 0;
k61 = k16;
k62 = k26;
k63 = k36;
k64 = k46;
k65 = k56;
k66 = [mm] 3*E2*I2*(l/2)^3/((l/2-l1)^3*l1^3);
[/mm]
%Steifigkeitsmatrix
K = [k11 k21 k31 k41 k51 k61; k12 k22 k32 k42 k52 k62; k13 k23 k33 k43 k53 k63; k14 k24 k34 k44 k54 k64
; k15 k25 k35 k45 k55 k65 ; k16 k26 k36 k46 k56 k66]
%Elemete der Dämpfungsmatrix
d11 = 1; % aus q1 = [1 0 0 0 0 0]
d12 = -1;
d13 = 0;
d14 = -1;
d15 = -1;
d16 = 0;
d21 = d12; % aus q2 = [0 1 0 0 0 0]
d22 = 1;
d23 = 0;
d24 = 1;
d25 = 1;
d26 = -1;
d31 = d13;
d32 = d23;
d33 = 1;
d34 = 1;
d35 = 0;
d36 = 0;
d41 = d14;
d42 = d24;
d43 = d34;
d44 = 1;
d45 = 1;
d46 = 0;
d51 = d15;
d52 = d25;
d53 = d35;
d54 = d45;
d55 = 1;
d56 = 0;
d61 = d16;
d62 = d26;
d63 = d36;
d64 = d46;
d65 = d56;
d66 = 1;
%Dämpfungsmatrix
D = [d11 d21 d31 d41 d51 d61; d12 d22 d32 d42 d52 d62; d13 d23 d33 d43 d53 d63; d14 d24 d34 d44 d54 d64
; d15 d25 d35 d45 d55 d65 ; d16 d26 d36 d46 d56 d66] ;
%Massenmatrix
M = [m1 0 0 0 0 0; 0 m2 0 0 0 0; 0 0 m3 0 0 0; 0 0 0 m4 0 0; 0 0 0 0 m4 0; 0 0 0 0 0 m6];
Mi = inv(M) ; %invertierte Massenmatrix
% Matrix F = Mi*K
Kr = Mi*K;
KrT = transpose (Kr);
% Matrix G = Mi*D
G = Mi*D;
GT = transpose (G);
% Shakerkraft
%Parameter der Kraft
amp = 9.1;
freq = 1
nullphase = pi/3
%Kraft (f) und Kraftvektor
f = amp*sin(2*pi*freq*t+nullphase)
F = [0;0;f;0;0;-f]
%Gleichungssystem übeführt in System erster ordnung
yp(1) = y(7)
yp(2) = y(8)
yp(3) = y(9)
yp(4) = y(10)
yp(5) = y(11)
yp(6) = y(12)
yp(7) = -((transpose (KrT*[1;0;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[1;0;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(8) = -((transpose (KrT*[0;1;0;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[0;1;0;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(9) = -((transpose (KrT*[0;0;1;0;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)])+transpose (Mi*[0;0;1;0;0;0])*F %-((transpose (GT*[0;0;1;0;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %+f %+shakerkraft
yp(10) = -((transpose (KrT*[0;0;0;1;0;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[0;0;0;1;0;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(11) = -((transpose (KrT*[0;0;0;0;1;0]))*[y(1);y(2);y(3);y(4);y(5);y(6)]) %-((transpose (GT*[0;0;0;0;1;0]))*[y(7);y(8);y(9);y(10);y(11);y(12)])
yp(12) = -((transpose (KrT*[0;0;0;0;0;1]))*[y(1);y(2);y(3);y(4);y(5);y(6)])+transpose (Mi*[0;0;0;0;0;1])*F %-((transpose (GT*[0;0;0;0;0;1]))*[y(7);y(8);y(9);y(10);y(11);y(12)]) %-f %-shakerkraft
end
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 12:20 Mi 22.08.2007 | | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|
