Verfasst am: 06.02.2013, 20:45
Titel: Modellberechnung mit Reibung
Hallo zusammen,
ich bräuchte wieder einmal Hilfe für die Lösung eines MATLAB-Problems und wäre für Ideen und Anregungen sehr dankbar!
Ich habe ein Feder-Masse-Modell, wobei sich die Masse auf einer Reibfläche bewegt. Die AB sind v(0)~=0 und s(0)=0. Mittels der Event-Funktion überprüfe ich bei jedem Stillstand der Masse (v=0), ob die Federkraft noch größer ist als die Reibkraft und sich die Masse somit weiterbewegt. Ich habe die DGL-Berechnung in eine Schleife eingebettet und sollte die Haftkraft größer als die Federkraft sein, ist die Bewegung und somit die Rechnung beendet. Hiermit sind wir bei meiner Frage:
Gibt es eine Möglichkeit, dass die Rechnung bis zum Ende der Simulationszeit weiterläuft, obwohl die Masse stillsteht?
Das Problem ist nämlich folgendes: Beende ich die Schleife (und damit die Rechnung) nicht bei v=0 und F_Feder<F_Haft, pendelt die Geschwindigkeit numerisch bedingt minimal um 0. Dadurch rechnet sich das Programm tot und natürlich springt die Beschleunigung stetig hin und her. Die Simulation kann somit nicht beendet werden. Kennt jemand eine Möglichkeit der Simulation für diesen Fall? Das Problem müsste eigentlich bei jeder Reibungssimulation auftreten. Das Weiterrechnen ist mir deshalb so wichtig, da dieses Reibungsmodell in ein größeres Modell implementiert werden soll, bei welchem die Rechnung nicht bei Stillstand der Batterie beendet werden darf.
Anbei findet Ihr die wichtigsten Teile meines Codes.
Ich hoffe sehr, dass Ihr mir weiterhelfen könnt.
while cont==1 [t y te ye ereig]=ode45(@SystemDGL,t_span,x0,options); % Aufruf der DGL-Function und Event-Function
tout=[tout; t]; % Ergebnisintervalle werden angehängt
yout=[yout; y];
lent=length(t); % Laenge des Vektors t (letzte Stelle gesucht) if(t(lent)>=t_e)% letzter Zeitwert größer Intervallende? disp(['Ende der Simulationszeit (t_e= ', num2str(t_e) ' s) erreicht']); % Ereignisausgabe im Command Window
cont=0; % -> Abbruch der while-Schleife else% Ereignis Grund fuer Abbruch? Ja, dann
telen=length(te); % Laenge des Vektors te (letzte Stelle gesucht)
F_c=abs(c_B*ye(telen,1)); % [N] aktuelle Federkraft
F_rH=(mue_H+mue_not)*m_B*g; % [N] (aktuelle) Reibkraft%
if ereig(telen)==1 && F_c<F_rH % Ereignis 1 eingetreten und Federkraft kleiner Reibungskraft disp(['Ereignis 1(v_B=0) und F_c<F_rH: Batteriestillstand nach t = ', num2str(te(telen)), ' s bei x = ', num2str(ye(telen)), ' m']); % Ereignisausgabe im Command Window
cont=0;
else disp(['Ereignis 1(v_B=0) und F_c>F_rH: Batteriestillstand nach t = ', num2str(te(telen)), ' s bei x = ', num2str(ye(telen)), ' m und bewegt sich weiter']); % Ereignisausgabe im Command Window
t_span=[te(telen) t_e]; % weiter mit Restintervall
x0=[ye(telen,1) ye(telen,2)]; % ... und neuen AB end end% Ende if/else-Bedingung end% Ende while-Schleife
Du bist der Experte, niemand wird mal eben verstehen was du da machst,
du musst die Frage etwas konkretisieren glaube ich.
Wie funktioniert die "Einbettung" in das übergeordnete System?
Die Schleife rechnet sich in jedem Fall tot wenn du sie nicht von
Hand abbrichst ob mit oder ohne Restgeschwindigkeit.
Und solange sie rechnet, kann auch nichts anderes rechnen
(Matlab hat nur einen Thread)
Und solange nichts anderes rechnet, wird sich auch nichts mehr ändern.
Soll doch noch was anderes rechnen, musst du es eben mit in die Schleife
hineinnehmen (nach der Federberechnung). Und bei v=0 setzt
du a einfach = 0 egal was die simulation ausrechnet
falls es noch Fragen gibt, bitte konkretisieren und in mögliche
Leser hineinversetzen!
Tut mir leid, wenn ich es zu ungenau beschrieben habe. Aber leider muss ich Dir auch widersprechen. Nehmen wir den gleichen Code wie oben und simulieren anstatt einem reibungsbehafteten Feder-Modell ein Feder-Dämpfer-Modell, so kommt die Bewegung / Masse auch zum Stillstand. Hier aber bleibt v=0 und die Rechnung läuft auch bis zum Ende der angegebenen Simulationszeit. Bei einem Reibungs-Modell springt v aber numerisch bedingt gering um 0 (je Iterationsschritt). Dies hat zur Folge, dass aufgrund μ*m*g*sign(v) die Beschleunigung ebenfalls zu jedem Zeitschritt hin und herspringt. Somit kommt die Bewegung / Masse nie zur Ruhe und ein Ende der Simulation ist nicht erreichbar. Ich möchte dem Programm irgendwie "sagen", dass er das v=0 annehmen soll, um die Rechnung beenden zu können.
Wie soll ich denn einfach a=0 setzen? Diese wird doch später aus DGL und den s- und v-Werten berechnet...
Und das v=0 taucht im Code nur in der Eventfunktion auf. Hierbei muss ja auch überprüft werden, wann die Masse zum Stillstand kommt, und nicht deren Beschleunigung. (Obwohl das in diesem Fall eh zum gleichen Zeitpunkt geschieht.)
Zwecks der Einbindung in das Gesamtmodell:
Hierbei meinte ich keine übergeordnete Funktion, sondern einfach ein größeres mechanisches System. Somit erweitern sich nur die DGL, die Berechnung bleibt aber fast diesselbe. Da dabei mehrere Massen beteiligt sind, ist es für die Berechnung des Systems wichtig, dass die Simulation nicht beim Stillstand einer Masse beendet wird.
Ich hoffe ich konnte das Problem besser erläutern. Falls nicht, versuche ich gern noch einmal mein Bestes. Zur Not erkläre ich auch nochmal das Gesamtmodell. Aber wichtig wäre erst einmal das reibungsbehaftete Feder-Modell zu lösen!
ok ich glaube ich fange an das Problem zu verstehen.
du löst die DLG ein Stück weit und gehst mit den Ergebnissen
als Anfangsbedingungen neu rein und löst wieder ein stück,
v wird nie ganz 0 und eine Mindestabbruchtoleranz kannst du nicht
festlegen, da sich das System ja an einem Umkehrpunkt befinden könnte.
D.h. es kommt definitiv auf die VORGESCHICHTE der Bewegung an,
also geht es um PERSISTENZ. Damit meine ich, es ist wichtig zu wissen
was in der längeren Vergangenheit passiert ist.
Du könntest
1. eine Variable einführen, die die GESCHICHTE der x_Position
aufzeichnet, und den Abstand ihrer Extremalwerte prüfen (Umkehrpunkte)
wenn diese Differenz unter einen bestimmten Wert sinkt, brichst du
die Schleife ab. Evt. musst du die ersten paar Schleifeniterationen
noch ignorieren, falls jeweilige Simulationsdauer sehr kurz ist.
Denk immer daran das Problem erst mal in Worten
zu formulieren sonst findest du nie eine programmatische Lösung.
Zb.:
> "Brich ab, wenn die Pendelbewegung winzig wird!"
> "Wie stelle ich fest wenn die Pendelbewegung winzig wird?"
> "Ich muss die Umkehrpunkte aufnehmen!"
2. du könntest das System so abändern, das dieses Verhalten nicht
mehr auftritt (vielleicht einen winzig kleinen Dämpfer oder so),
das sind aber fachliche Dinge die nur du entscheiden kannst
du hast das System und das Problem jetzt genau verstanden!
Deine Idee mit der Überprüfung der letzten Werte ist mir auch schon eingefallen. Aber ich will ja dann die Schleife NICHT abbrechen, sondern irgendwie realisieren, dass er mit v=0 trotzdem noch bis zum Ende der Simulationzeit rechnet!
naja dann eben nicht abbrechen, sondern anders reagieren,
aber du musst erstmal feststellen DASS dieser ungeliebte Fall
überhaupt eingetreten ist, Matlab sagt es dir nicht.
dazu kannst du die vorgeschlagenene Methode verwenden.
Ohne das Gesamtsystem zu kennen werde ich niemals verstehen was
überhaupt gemeint ist mit "weiterrechnen". Kann ein weiteres Systemglied
dein Feder-Reibsystem wieder neu anregen nachdem es schon
still stand (bzw Mikropendelt)?
Und wenn ja, dann musst du IN der Differentialgleichung verzweigen:
Stellst du ein Mikropendeln fest dann schaltest du in den Pauschalzustand
v = 0. Übersteigt die Federkraft später wieder die Reibkraft schaltest du wieder in den normalzustand.
Also eine Differentialgleichung die sich dynamisch anpasst was mit
einer if-Verzweigung leicht machbar sein sollte
schon einmal Danke, dass du dich dem Problem so annimmst!
Zwecks des Weiterrechnens: Ich hab in das System beispielhaft einen Kraftimpuls bei 3s für 0,001s mit der Eventfunktion eingebaut. Funktioniert auch. Somit kann ich den gewünschten Effekt überprüfen.
Deine Idee mit dem Ab- und Zuschalten ist sehr gut und das hatte ich auch schon versucht. Nur weiß ich nicht genau wo und wie ich das machen soll. Bisher habe ich in der DGL-Funktion folgendes gemacht:
Die Konstante p stellt sich bei v=0 und F_Feder<F_Haft automatisch auf 0, damit die DGL mit v=0 weiterberechnet wird. Sobald F_Feder>F_Haft, stellt sich p=1 und die normale DGL wird weiterberechnet. Das klappt auch für die s- und v-Werte sehr gut, diese werden "eingefroren". Das Problem ist, dass der a-Wert auch "eingefroren" wird, aber leider nicht bei 0, sondern bei a~=0.
Um das zu verdeutlichen, siehst du im Anhang den Ergebnisplot.
Hast Du evtl. eine bessere Lösung? Wie stellst Du Dir die if-Verzweigung vor? Ich hatte schon auf eine saubere Lösung als die meine oben gehofft.
, spielen hier s_B und v_B die entscheidende Rolle.
s_B) Da die Masse auch bei s_B~=0 stehen bleiben kann, muss ich s_B mit p*s_B auf 0 setzen, da hast du Recht.
v_B) Hier das Problem. Wenn ich in der DGL-Funktion v_B auf 0 setze, dann müsste doch im Hauptprogramm bei y bzw yout auch 0 "ankommen". Das tut es aber nicht, eher ist v_B~=0 -> sign(v_B)~=0 -> a_B~=0. Ist das richtig? Habe ich hier ein programmtechnischen Denkfehler?
Sicherlich kann ich auch p*v_B auch in die Beschleunigungsberechnung integrieren und auch im Gesamtmodell ausprobieren.
Trotzdem würde ich gern verstehen, warum die v-Wert zwischen DGL-Funktion und Hauptfunktion unterschiedlich sind.
Das werden saubere Nullen. Ich gehe davon aus, dass du p auf 0
setzt wenn v eine sehr kleine Geschwindigkeit unterschreitet, weil
du genau 0 vielleicht nicht erreichst. Diese winzige Geschwindigkeit
geht aber als Anfangsbedingung wieder in den nächsten Rechenzyklus
und innerhalb der DGL hast du keinen Einfluss darauf, weil die
DGL ja ausrechnen soll was sich an diesem Zustand ÄNDERT
Und bei p = 0 ändert sich...... nichts!
Die winzige Geschwindigkeit der Anfangsbedingung bleibt schlicht bestehen.
Die Lösung ist dann natürlich, bei p = 0 auch die neuen
Anfangsbedingungen auf 0 zu setzten oder eben den "falschen"
Output wieder mit p zu multiplizieren (wie du vorgeschlagen hast)
Das wäre meine Theorie. Gib mal bitte ein Feedback
deine Idee ist sehr gut und vollkommen plausibel. Um dies zu überprüfen habe ich in der Schleife für alle neuen Geschwindigkeitsanfangswerte "p*" eingebaut.
> x0=[ye(telen,1) p*ye(telen,2)];
Dann erhalte ich die Fehlermeldung:
??? Error using ==> horzcat
CAT arguments dimensions are not consistent.
Error in ==> ode45 at 425
ieout = [ieout, ie];
Error in ==>
[t y te ye ereig]=ode45(@SystemDGL,t_span,x0,options);
Ich schlußfolgere, dass MATLAB mit zwei zum gleichen Zeitschritt auftretenden Events nicht umgehen kann, hier v=0 und t=t_Kraftimpuls. Ich habe auch keine Lösung gefunden, die ein solches Eventproblem angeht.
Testweise habe ich den Kraftimpuls deaktiviert. Der gewünschte Effekt ist sichtbar. v ist wirklich 0 und somit die Beschleunigung auch exakt 0 (ohne "p*" in der Beschleunigungsformel) und die Simulation läuft bis zum vorgegebene Ende.
Ich werde nun beide Methoden
1) "p*" in der Beschleunigungsformel
2) "p*" in den AB's
im Gesamtmodell ausprobieren.
Sollte Dir eine Lösung zu der Event-Problematik einfallen, würde ich mich über ein Feedback freuen.
hm naja, die Event-Konflikt-Theorie kann ich mir irgendwie nicht
so richtig vorstellen. Ich glaube da überschätzt du Matlab ^^.
ML gibt ja blos stumpf Argumente weiter, was ein "Event" ist und
was nicht ergibt sich erst aus der physikalischen Interpretation der Daten.
1. Vorschlag:
Setzte Breakpoints nach x0 = [ ... p* ...] und zwar im Code mit
"p*..." und ohne "p*..." und schaue dir x0 an. Haben die das gleiche Format?
2. Vorschlag:
poste mal die beiden codes. Ich hab seit langer Zeit nicht mehr den
aktuellen Stand und muss viel spekulieren.
die Fehlermeldung besagt, dass die Dimensionen von ie, also dem Ereignisvektor, überschritten werden. Das kann ja nur sein, wenn mehr als ein Event in einem Zeitschritt auftritt und die Infos somit nicht im ie-Vektor abgespeichert werden können. Dies wird auch dadurch bestätigt, dass wenn ich die Kraftimpuls-Events rausnehmen, das Programm ohne Fehlermeldung durchläuft. Außerdem tritt der Error direkt beim Kraftimpuls_Beginn-Event auf.
x0 ohne p) x0=[... ~=0]
x0 mit p) x0=[... =0]
-> deine Idee mit den AB's war somit vollkommen richtig.
Der derzeitige Code mit Krafimpuls und p*AB->der Code mit Errror:
zunächst mal muss ich eingestehen, dass ich mit den ode-events
nicht 100% vertraut bin
Aber nach dem was ich jetzt gelesen habe, gibt es keine konkurrierenden
Events. Es gibt blos das eine "Nulldurchgang"-Event. Und da dann
die 3 Arten: direction = Von oben, Von unten, Egal
Die genaue Fehlermeldung lautet ja: horzcat error
Dem, und aus einem Beispiel in der Hilfe folgend würde ich mal
versuchen: value, direction und isterminal als SPALTENVEKTOREN
anzulegen
value = [ v_b ; wasweisich ]
Ich kann es leider nicht testen da dein geposteter Code so nicht
lauffähig ist (es fehlen odeset, p, die Anfangsbedingungen usw).
Ein lauffähiger Code ist immer um ein vielfaches hilfreicher.
Wenn das nicht hilft, nochmal melden, das interessiert mich jetzt auch.
also ich habe keine Ahnung warum, aber es hat funktionert. Ich habe [value,isterminal,direction] anstatt als Zeilenvektoren nun als Spaltenvektoren abgespeichert und übergeben und das Programm läuft trotz zeitgleiches Event-Auftreten durch.
Danke für die Idee!!!
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS