Anwendung unterschiedlicher ode-Solver am realen Beispiel

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

JayvH

Forum-Anfänger


	Beiträge: 36

	Anmeldedatum: 19.03.11

	Wohnort: Hannover

	Version: ---

Verfasst am: 03.06.2012, 15:17 Titel: Anwendung unterschiedlicher ode-Solver am realen Beispiel

Hallo gomatlab-Forum,

ich versuche derzeit die Wellenverlagerungsbahn eines Gleitlagers zu bestimmen. Ich glaube ich habe es mittlerweile auch geschafft, bin aber bei manchen etwas extremen Geometrien dazu gezwungen anderen ode-Solver zu verwenden bzw. bekomme den Standard-Solver ode45 nur mit einigen Verrenkungen dazu, zu tun was ich möchte. Ich überblicke leider nicht, warum gerade der Solver im Gegensatz zu den anderen nicht so leicht funktioniert. In der Literatur bzw. diversen Internetseiten zu dem Thema heißt es ja öfter mal etwas zu steifen und nicht steifen Differentialgleichungen, wo ich einfach nicht weißt, ob das auf mein System zutrifft. Es ist ein bißchen Fleißarbeit das alles einzugeben, aber ich hoffe, damit wird es etwas leichter mein Problem zu beschreiben und vielleicht kann mir auch jemand etwas auf die Sprünge helfen, damit ich mich durch MATLAB nicht mit Trial & Error kämpfe, sondern auch noch etwas dabei lerne.

Diese wird bestimmt durch die relative Lage zur Lagermitte $\epsilon$ sowie den Winkel zur Senkrechten $\delta$ .

Die zeitlichen Veränderungen sind definiert durch folgende Differentialgleichungen:
$\dot{\epsilon} = \frac{F \cdot \psi^2}{b \cdot d \cdot \eta \cdot So_{V}} \cdot \left[ \cos (\delta - \gamma) - \frac{sin | \delta - \gamma|} {\tan \beta} \right]$
$\dot{\delta} = \frac{\omega_{W} + \omega_{W}}{2} - \frac{F \cdot \psi^2}{b \cdot d \cdot \eta \cdot So_{D}} \frac{sin (\delta - \gamma)} {2 \sin \beta}$

Die Parameter $So_{V}$ , $So_{D}$ und $\beta$ werden durch Näherungsfunktionen bestimmt. Sie sind alle abhängig von $\epsilon$ sowie der Lagerbreite $b$ und dem Durchmesser $d$ und lauten im einzelnen:
$So_{D} = \left(\frac{a}{b} \right)^2 \cdot \frac{\epsilon}{2 \left( 1 - \epsilon^2 \right)^2} \sqrt{\pi^2 \left(1 - \epsilon^2 \right) + 16 \epsilon^2} \cdot \frac{a_1 \cdot \left(\epsilon - 1 \right)}{a_2 + \epsilon}$
$\beta = \arctan \left( \frac{\pi \sqrt{1 - \epsilon^2} }{2 \epsilon} \right) \cdot \left( a_3 + a_4 \epsilon + a_5 \epsilon^2 + a_6 \epsilon^3 + a_7 \epsilon^4 \right)$
$So_{V} = 4 \left(\frac{b}{d}\right)^{2} \left(1 - \epsilon^2\right)^{-\frac{5}{2} } \left[ \left( \frac{\pi}{2} - \frac{1}{2} \arccos \epsilon \right) \left(1 + 2\epsilon^2 \right) + \frac{3}{2} \epsilon \sqrt{1 - \epsilon^2} \right] \cdot \frac{a_{8} \left(1-\epsilon\right)}{-a_{9} - \epsilon}$ , wo $\epsilon$ negativ einzusetzen ist, wenn die Verdrängung negativ ist. Das alles habe ich in eine MATLAB-Funktion integriert, die so aussieht:

Code:

function dot_y = DGL_System_Lang(t, y, F, gamma, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, b, d, psi, eta, omega_W, omega_L, Verdraengung_negativ)
epsilon = y(1);
delta = y(2);
dot_y = zeros(2,1);
dot_y(1) = (F * psi^2)/(b * d * eta * fun_So_V(epsilon, a_8, a_9, b, d, Verdraengung_negativ)) * (cosd(delta - gamma) - (sind(abs(delta - gamma)))/(tand(fun_beta(epsilon, a_3, a_4, a_5, a_6, a_7))));
dot_y(2) = radtodeg((omega_W + omega_L)/2 - (F * psi^2)/(b * d * eta * fun_So_D(epsilon, a_1, a_2, b, d)) * (sind(delta - gamma))/(2 * sind(fun_beta(epsilon, a_3, a_4, a_5, a_6, a_7))));

% Ermittlung Sommerfeldzahl für Drehung
function So_D = fun_So_D(epsilon, a_1, a_2, b, d)
So_D = (b/d)^2 * epsilon/(2 * (1 - epsilon^2)^2) * sqrt(pi^2 * (1 - epsilon^2) + 16 * epsilon^2) * (a_1 * (epsilon - 1))/(a_2 + epsilon);

% Ermittlung Sommerfeldzahl für Verdrängung
function So_V = fun_So_V(epsilon, a_8, a_9, b, d, Verdraengung_negativ)
if strcmp(Verdraengung_negativ, 'false') % Fallunterscheidung
So_V = 4 * (b/d)^2 * (1 - epsilon^2)^(-5/2) * ((pi/2 - 1/2 * acos(epsilon)) * (1 + 2 * epsilon^2) + 3/2 * epsilon * sqrt (1 - epsilon^2)) * (a_8 * (1 - epsilon))/(-a_9 - epsilon);
else
epsilon = -epsilon;
So_V = 4 * (b/d)^2 * (1 - epsilon^2)^(-5/2) * ((pi/2 - 1/2 * acos(epsilon)) * (1 + 2 * epsilon^2) + 3/2 * epsilon * sqrt (1 - epsilon^2)) * (a_8 * (1 - epsilon))/(-a_9 - epsilon);
end

%Ermittlung beta
function beta = fun_beta(epsilon, a_3, a_4, a_5, a_6, a_7)
beta = atand((pi * sqrt(1 - epsilon^2))/(2 * epsilon)) * (a_3 + a_4 * epsilon + a_5 * epsilon^2 + a_6 * epsilon^3 + a_7 * epsilon^4);

Funktion ohne Link?

Der Aufruf erfolgt dann in meinem Hauptskript, wo ich auch die Eingaben für die anderen Parameter vorgebe über folgenden Aufruf:

Code:

if epsilon(Zaehler) > 0.95
Genauigkeit = 'hoch';
else
Genauigkeit = 'normal';
end
switch Genauigkeit
case 'hoch'
options = odeset('maxstep',1e-7);
% options = odeset('AbsTol', 1e-48, 'RelTol', 1e-24);
[t, y] = ode45(@(t, y) DGL_System_Lang_Interpolation(t, y, Kraftfunktion, Winkelfunktion, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, b, d, psi, eta, Delta_t, omega_W(Index), omega_L(Index), Zaehler), [Zeit(Zaehler) Zeit(Zaehler+1)], [epsilon(Zaehler) delta(Zaehler)], options);
case 'normal'
[t, y] = ode45(@(t, y) DGL_System_Lang_Interpolation(t, y, Kraftfunktion, Winkelfunktion, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, b, d, psi, eta, Delta_t, omega_W(Index), omega_L(Index), Zaehler), [Zeit(Zaehler) Zeit(Zaehler+1)], [epsilon(Zaehler) delta(Zaehler)]);
end

Funktion ohne Link?

Wie man sieht brauche ich hier für den ode45-Aufruf eine Fallunterscheidung. Falls das Verhältnis von $\frac{b}{d}$ sehr klein wird (in meinem Fall ein Achtel), dann wird $So_V$ sehr klein und $\dot{\epsilon}$ sehr groß. $\epsilon$ kann dadursch schnell größer als 1 werden. Was rein physikalisch nicht möglich ist, da das Lager dann durch die Wand geht und hier in der Formel für $\beta$ bestraft wird, da der Wurzelausdruck komplex wird und der arctan dann nicht mehr funktioniert. Durch eine Verkleinerung der stepsize kann ich das hinbekommen, aber der Rechenvorgang dauert dann sehr lange. Wie man sieht habe ich mich auch schon an der Einstellung über die Toleranzen versucht, was ja das probatere Mittel sein soll. Leider habe ich das aber gar nicht hinbekommen. Die Verwendung eines jeden anderen Solvers vermeidet das Problem.

Frage(n):
Handelt es sich hierbei um ein steifes System, weshalb die anderen Solver besser sind? Wie könnte man das Problem vernünftig über die Einstellung der Fehlertoleranzen lösen ohne die Brachiallösung der Maxstepsize nehmen zu müssen?

Wenn es hilft, dann könnte auch mal ein vollständiges Beispiel mit Kräften und Geometriedaten hochladen.

Wäre nett, wenn mir jemand ein paar Hinweise zu besseren Verständnis geben könnte.

JayvH

Themenstarter

Forum-Anfänger


	Beiträge: 36

	Anmeldedatum: 19.03.11

	Wohnort: Hannover

	Version: ---

Verfasst am: 18.06.2012, 19:35 Titel:

Bis jetzt leider keine Antwort. Kann sein, dass das Beispiel eher abschreckt als zu Diskussionen anregt. Von daher vielleicht mal losgelöst davon:

Gehe ich recht in der Annahme, dass eine steife DGL bedeutet, dass die Änderung der Zustandsgröße in Abhängigkeit eines Parameters extrem stark ist?


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	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

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.