Grundlagen: Simulation - Solver

Autor: Dr.-Ing. Alfred Kamusella

Simulationen mit numerischen Modellen beruhen im Wesentlichen auf der Berechnung von Modell-Ergebnissen mit Hilfe eines Gleichungslösers (Solver). In diesem Grundlagen-Kapitel liegt der Schwerpunkt auf der Erläuterung der prinzipiellen Wirkungsweise von Solvern zur Simulation dynamischer Systeme im Zeitbereich. Dabei wird beispielhaft Bezug genommen auf die Begriffe und Konzepte, wie sie im Simulationssystem SimulationX verwendet werden.

Numerische Simulation zeitkontinuierlicher Systeme

Kenngrößen eines Simulationslaufes (Zeitachse)

Aus Sicht des Modells läuft die Simulationszeit time bzw. t in diskreten Schritten dt von tStart bis tStop (meist vorwärts):
Innerhalb des aktuellen Zeitintervalls dt werden in Abhängigkeit vom verwendeten Lösungsverfahren mehrere Stützstellen berechnet:
- Jede Stützstelle entspricht einer einmaligen Modellberechnung. Demzufolge sind die meisten Modelldurchrechnungen Hilfsrechnungen (mit "unexakten" Werten).
- Nach Durchführung aller Hilfsrechnungen erfolgt einmalig für das Ende des aktuellen Zeitintervalls dt eine Modellberechnung mit den "richtigen" Modellwerten.
- Innerhalb des Modell-Algorithmus kann man auf die letzten "richtigen" Modellwerte zugreifen (in SimulationX mit der last-Funktion).
Bei numerischen Verfahren mit Schrittweitenregelung erfolgt eine automatische Anpassung der Schrittweite dt an die Eigenschaften des Modells zum berechneten Zeitpunkt ("Wie schnell können sich die physikalischen Zustandsgrößen ändern?"):
- Je empfindlicher Zustandsgrößen des Modells auf andere Modellgrößen reagieren können, desto kleiner wird dt gewählt.
- Die zulässige Rechenschritte dt kann man begrenzen (dtMin soll "ewiges" Rechnen vermeiden / dtMax soll eckige Signalverläufe vermeiden).
Die Berechnung mit der sich kontinuierlich anpassenden Schrittweite dt wird gestört durch Ereignis-Zeitpunkte, welche "exakt" angesteuert und berechnet werden müssen. Die Wirkungsweise der erforderliche Ereignisbehandlung wird in einem gesonderten Abschnitt erläutert.

Abschnitte eines Simulationslaufes

Ein Simulationslauf im Zeitintervall tStart ≤ time ≤ tStop besteht aus den Abschnitten INIT, DYNAMIC und FINISH. Im Folgenden wird die Aufgabe dieser Abschnitte zum einen aus numerischer Sicht erläutert. Zum anderen erfolgen Hinweise, welche Aufgaben die zugehörigen Modellabschnitte für die Nutzung von Modellen im Konstruktionsprozess übernehmen:

Grundlagen Simulation - Modellberechnung - zeitkontinuierliche Systeme - INIT-Abschnitt - beispiel ode.gif

INIT ( time = tStart )
dient der Ermittlung konsistenter Anfangswerte für alle Zustandsgrößen (Bilder aus der SimulationX-Hilfe):

Die Anfangswerte charakterisieren den (energetischen) Zustand des Modells zu Beginn der Simulation.
Anfangswerte gewöhnlicher DGL sind immer konsistent (ODE: Ordinary Differential Equation). Der Solver kann daraus eindeutig die Werte der höchsten Ableitungen berechnen (im Beispiel Beschleunigung a(tStart)=-1100 m/s²).
Modelle mit Zwangsbedingungen (z.B. für Getriebe) führen zu Algebro-Differenzialgleichungen (DAE: Differential Algebraic Equation). Hier können einige Zustandsgrößen nicht unabhängig voneinander gewählt werden:
Man muss sich entscheiden, welche Anfangswerte man definiert vorgibt (fixiert). Die nicht fixierten Anfangswerte versucht der Solver, aus den im Modell in Form von Gleichungen beschriebenen Zwangsbedingungen zu berechnen. Dies gelingt insbesondere bei stark nichtlinearen oder nicht nichteindeutigen Zusammenhängen nicht immer!
Man sollte möglichst viele gültige Anfangswerte vorgeben und fixieren! Nur dort, wo die analytische Berücksichtigung der Abhängigkeiten praktisch unmöglich ist, überlässt man dem Solver die Anfangswertbestimmung!

Idealisierte Netzwerkelemente besitzen konzentrierte Parameter (Punktmasse, Elastizität, Widerstand usw.):

Über die Dimensionierungsgleichungen der Elemente ist eine Transformation konstruktiver Kenngrößen (Geometrie und Stoffeigenschaften) in die erforderliche konzentrierten Parameter möglich.
Für die zeitunabhängigen konzentrierten Parameter sollte diese Dimensionierungsgleichungen in den INIT-Abschnitt der Modell implementiert werden. Es entstehen dadurch Modell-Elemente, welche mit konstruktiven (stofflich-geometrischen) Parametern gespeist werden (Abmessung, Dichte, E-Modul usw.).

DYNAMIC ( tStart ≤ time ≤ tStop )
Es wirkt die Überführungsfunktion des dynamischen Systems (DGL=dynamisches Modell):

Die numerische Lösung des DGL-Systems berechnet die zeitliche Entwicklung aller Zustandsgrößen in Form von Signalverläufen.
Alle Modellgrößen, deren Werte von den Zustandsgrößen abhängig sind, werden ebenfalls als Signalverläufe berechnet.
Die mit dem Modell berechneten Signalverläufe repräsentieren das zeitliche Verhalten des abgebildeten Originals.

FINISH ( time = tStop )
ermöglicht die abschließende Berechnung integraler Werte (z.B. mittlere Verlustleistung, Statistik, Fourieranalyse):

Im FINISH-Abschnitt des Modells sollte die automatisierte Auswertung der im DYNAMIC-Abschnitt berechneten Signalverläufe implementiert werden.
Wichtig für die weitere Verwertung der Simulationsergebnisse im Konstruktionsprozess sind Bewertungsgrößen, welche die Qualität des Modellverhaltens repräsentieren (z.B. Wirkungsgrad oder Kosten aus Material- und Energieverbrauch). Möglich wäre hier auch eine "ästhetische" Bewertung im Sinne der Eleganz des Bewegungsablaufs in einem Antrieb.

Ableitungs- und Integralform von DGL

Die Überführungsfunktion eines dynamischen Modells stellt ein DGL-System dar. Zur Lösung dieses DGL-Systems verwenden die Solver Verfahren der numerischen Integration:

Zweck der Verfahren zur numerischen Integration ist die näherungsweise Berechnung der Fläche (bestimmtes Integral) unter einem Funktionsverlauf (Ableitung).
Es existieren zwei grundsätzliche Möglichkeiten, um ein DGL-System zu notieren (hier am Beispiel des Pendels):
1. Ableitungsform:
2. Integralform:
Die Integralform enthält explizit die Zustandsgrößen als Ergebnis einer numerischen Integration. Insofern ist diese Form für den technischen Anwender meist anschaulicher als die Ableitungsform.
Unabhängig von der im Modell gewählten Notation erfolgt durch das Simulationssystem für den Anwender unsichtbar eine Transformation in die für den Solver geeignete Form eines DGL-Systems.

Wechselwirkung von Zustandsgroesse Y und Ableitung YP

Grundlagen Simulation - Modellberechnung - zeitkontinuierliche Systeme - flaechenintegral.gif

Y(t) als Resultat einer Integral-Bildung kann geometrisch als Fläche unter dem Funktionsverlauf der zugehörigen Ableitung YP(t) interpretiert werden:

Grundlagen Simulation - Modellberechnung - zeitkontinuierliche Systeme - Wechselwirkung Y und YP - integral.gif

Y(t) repräsentiert in einem dynamischen Modell die Füllung eines Speicher-Elements mit Energie, Stoff bzw. Information (=Zustandsgröße). YP(t) ist dabei die Änderungsgeschwindigkeit dY/dt für die Füllung des Speicherelements.
Im Allgemeinen ist YP(t)=f(Y(t)), d.h. die Geschwindigkeit des Energie- /Stoff- /Info- Flusses zwischen den Speichern ist abhängig von den Potentialdifferenzen und den Übertragungselementen zwischen den Speichern.
Da YP(t) somit nicht explizit gegeben ist, kann Y(t) nicht hinreichend genau durch numerische Verfahren der Flächenberechnung ermittelt werden (z.B.: Trapez- /Tangentenformel oder Simpsonsche Regel).
Die konkreten zeitlichen Verläufe von Y(t) und YP(t) ergeben sich erst durch die zeitliche Wechselwirkung zwischen den Energiespeichern!

Zentraler Integrationskern (Solver)

Es existieren verschiedene Verfahren zur Lösung (=Simulation) von DGL mittels numerischer Integration. Diese Verfahren sind meist in einem "zentralen Integrationskern" = "Solver" implementiert:

Die Integral- bzw. Differential-Funktion ist die Grundlage für die Verbindung zwischen Modell und Solver:

:
a  := F/(m);
v  := integral (a,v0);
x  := integral (v,x0);

Es sieht auf den ersten Blick so aus, als würde über diese Funktionsaufrufe der zentrale Integrationskern vom Modell aus aufgerufen, jedoch verhält es sich umgekehrt!
Das Modell ist aus Sicht der numerischen Integration eine Black-Box, welche die Ableitungen der Zustandsgrößen berechnet:

Prinzip der numerischen Integration

Das DGL-System wird durch die numerische Integration zum "Leben" erweckt. Es entwickelt eine Eigendynamik auf der Basis seiner Elemente (Speicher, Übertragungsglieder und Wandler) und ihrer Wechselwirkung. Der Solver ist in der Lage, über einen begrenzten zeitlichen Prognosehorizont den Zustand des Systems hinreichend genau vorauszusagen:

Ausgehend von den Werten der Zustandsgrößen und Eingangssignale zum aktuellen Zeitpunkt t_n sowie den Modellparametern werden mit dem Modell die zeitlichen Ableitungen YP(t) der Zustandsgrößen Y(t) ermittelt.
Mit dieser "Abschätzung" der zeitlichen Änderung der Zustandsgrößen prognostiziert der Solver über das Zeitintervall dt die Werte der Zustandsgrößen für den Zeitpunkt t_n+dt.
Im obigen Bild wird dies am Beispiel des EULER-Verfahrens verdeutlicht. Bei diesem Verfahren wird vereinfachend angenommen, dass die Ableitung der Zustandsgrößen im Zeitintervall dt konstant auf dem Wert YP(t) bleibt.
Nach der Art und Weise, wie die Prognose der Zustandsgrößen für den nächsten Integrationsschritt dt erfolgt, kann man die Verfahren der numerischen Integration grob in grundsätzliche Klassen einteilen. Diese Klassifizierung der Integrationsverfahren soll in den folgenden Abschnitten kurz vorgestellt werden.

Explizite und implizite Verfahren

Zum aktuellen Zeitpunkt t_n können mit dem dynamischen Modell die aktuellen zeitlichen Ableitungen der Zustandsgrößen berechnet werden:

Der effektive Wert der Ableitung der Zustandsgrößen gemittelt über das gesamte Intervall dt wird unabhängig von der Einteilung in explizit/implizit durch Hinzunahme weiterer Informationen "geschätzt". Die Prinzipien dieser Abschätzung werden in den folgenden Abschnitten noch erläutert.

Explizite Verfahren betrachten den Schätzwert für die effektive Ableitung YP als hinreichend genau. Grafisch stellt sich die Prognose einer Zustandsgröße für den Zeitpunkt t_n+dt als Anstieg der Tangente YP an die Zustandsgröße YP(t_n) im Intervall dt dar:

Der prognostizierte Wert der Zustandsgrößen wird dann ebenfalls als hinreichend angenommen und ist Ausgangspunkt für die Berechnung des nächsten Zeitschrittes dt.

Implizite Verfahren bewegen sich im Zeitbereich nicht so geradlinig vorwärts wie explizite Verfahren. Die Berechnung jedes Zeitschrittes dt ist durch eine Iterationsschleife gekennzeichnet:

Unter Berücksichtigung des berechneten Anstiegs YP(t_n) wird ein Anfangswert für die Zustandsgröße am Ende des Zeitschrittes dt Y(t_n+dt) abgeschätzt (z.B. durch Extrapolation des bisherigen Signalverlauf Y(t)).
Der mit dem aktuellen Wert von Y(t_n+dt) mittels des dynamischen Modells berechnete Anstieg YP(t_n+dt) wird als Tangente an Y(t_n+dt) angelegt und rückwärts zum aktuellen Zeitpunkt t_n verfolgt.
Ziel der anschließenden Iteration ist die Verbesserung des Schätzwertes für Y(t_n+dt) solange, bis die rückwärts verfolgte Tangente YP den Wert von Y(t_n) hinreichend genau trifft.
Ist dieses Ziel erreicht, werden die aktuellen Werte der Zustandsgrößen Y(t_n+dt) als hinreichend genau angenommen und das Integrationsverfahren beginnt mit der Berechnung des nächsten Zeitschrittes dt.

Einschritt- und Mehrschritt-Verfahren

Unabhängig von der Einteilung in explizite und implizite Verfahren existiert eine Unterteilung in die sogenannten Einschritt- und Mehrschritt-Verfahren.

Einschritt-Verfahren

"Betrachten" für die Prognose des nächsten Systemzustandes (Zeitpunkt=t_n+dt) ausgehend von Y(t_n) nur den aktuell auszuführenden Zeitschritt dt.
Je nach Genauigkeitsordnung des konkreten Verfahrens wird eine Anzahl von Systemabtastungen innerhalb dieses Zeitschrittes durchgeführt (=Modellrechnungen).
Nach der erfolgreichen Prognose für Y(t_n+dt) rückt das vom Verfahren betrachtete Zeitfenster um dt weiter.
Da nur Informationen aus dem aktuellen Simulationsschritt verwendet werden, sind starke Änderungen der dynamischen Beziehungen relativ unproblematisch (Schrittweitenregelung).

Mehrschritt-Verfahren

Durch Bestimmung eines Interpolationspolynoms, das durch zurückliegende Zeit-Stützpunkte der Lösungsfunktion gelegt wird, kann deren weiterer Verlauf prognostiziert werden.
Es wird zwischen reinen Mehrschrittverfahren (Adams-Bashforth-Verfahren) und den Prädiktor-Korrektor-Verfahren unterschieden.
Bei Prädiktor-Korrektor-Verfahren wird der extrapolierte Wert der Lösungsfunktion in mehreren Iterationen korrigiert (implizites Verfahren).
Bei stetigen Modellen benötigen Mehrschritt-Verfahren weniger Abtastungen als Einschritt-Verfahren und sind somit schneller.
Bei Unstetigkeiten geht dieser Vorteil durch die ständigen Wiederanlaufphasen schnell verloren (erste Schritte mit Einschritt-Verfahren!).

Beispiele fuer Integrationsverfahren

Im Folgenden soll exemplarisch die Wirkungsweise wichtiger Verfahren der numerischen Integration vorgestellt werden. Aus der unterschiedlichen Wirkungsweise resultieren die unterschiedlichen Eigenschaften der Integrationsverfahren. Dieses Wissen ist nützlich, wenn man für ein konkretes Modell ein geeignetes Verfahren der numerischen Integration auswählen und konfigurieren muss.

Explizite Einschritt-Verfahren

Das EULER-Verfahren ist das einfachste Verfahren:

Man nimmt an, dass die Änderungsgeschwindigkeit YP der Zustandsgrößen Y über die Zeitspanne dt hinweg konstant bleibt:
Infolge der Aufsummierung der Prognosefehler wird das EULER-Verfahren sehr schnell instabil!

Beim HEUN-Verfahren wird ein Korrekturschritt angefügt:

Mit der nach EULER prognostizierten Zustandsgröße Y₁(t+dt) berechnet man die Ableitung YP(Y₁(t+dt)).
Dann nimmt man den Mittelwert beider Ableitungen als genauere Annahme für den Anstieg YP(t):
Der Berechnungsaufwand für einen Zeitschritt dt verdoppelt sich im Vergleich zum EULER-Verfahren.
Infolge der verbesserten "Abschätzung" der effektiven Ableitung im Zeitintervall dt ist jedoch der Prognosefehler um Größenordnungen geringer als beim EULER-Verfahren und das Verfahren arbeitet auch mit größeren Schrittweiten dt wesentlich stabiler als das EULER-Verfahren.

Verallgemeinerung für explizite Runge-Kutta-Verfahren :

Das EULER- und HEUN-Verfahren sind Runge-Kutta-Verfahren 1. bzw. 2. Ordnung.
Je nach Ordnung des konkreten Runge-Kutta-Verfahrens wird der Funktionswert YP(t) auf der Basis mehrerer Stützstellen (Modellrechnungen) im betrachteten Zeitschritt "geschätzt":
Die Zahl der erforderlichen Stützstellen steigt mit der Ordnung des Verfahrens. Mit der Ordnung des Verfahrens steigt auch die Genauigkeit der Berechnung.
Durch Ermittlung von YP(t) mittels zweier unterschiedlicher Genauigkeitsordnungen kann eine Fehlerabschätzung und Schrittweitenregelung erfolgen. Man spricht hierbei von eingebetteten Verfahren, weil für das Verfahren der niedrigeren Ordnung eine Teilmenge der Stützstellen der höheren Ordnung verwendet werden.

Implizite Mehrschritt-Verfahren

Bei den BDF-Verfahren (Backward-Differential Formulas) handelt es sich um implizite Verfahren mit einstellbarer Ordnung und automatischer Schrittweitensteuerung:
Auf Grundlage des bereits berechneten Verlaufs der Zustandsgröße Y(t) bis zum Zeitpunkt t erfolgt eine Startwert-Ermittlung Y_o(t+dt) als Ausgangspunkt für die "exakte" Lösung Y(t+dt) je nach Ordnung des Verfahrens:
Danach erfolgt eine iterative Lösung des nichtlinearen Gleichungssystems zur Berechnung der Zustandsgrößen Y(t+dt) mit Abbruch bei der geforderten Genauigkeit.

Die BDF-Verfahren sind sogenannte Prädiktor-Korrektor-Verfahren:

Für jede Zustandsgröße wird bei einer Verfahrensordnung k aus k+1 bereits berechneten Werten der nächste aktuelle Wert als Anfangswert Y_o(t+dt) extrapoliert (Prädiktor = Prognose durch explizites Mehrschrittverfahren).
Der Wert Y(t+dt)wird iterativ durch ein implizites Mehrschrittverfahren gleicher Ordnung so lange verbessert, bis bestimmte Konvergenzkriterien erfüllt sind (corrector = Verbesserer):
In jeder Korrektor-Iteration muss ein lineares Gleichungssystem gelöst werden:

res = J · (Y_i-1 - Y_i)

mit res = Residuen / J = Jacobimatrix / Y = Zustandsgrößen

Das Residuum widerspiegelt den momentanen Fehler in den Differentialgleichungen.
Nach jeder Iteration wird die zu Y_i gehörende Ableitung mittels Modellberechnung ermittelt und "rückwärts" überprüft, ob damit Y_i-1 erreicht wird (implizites Verfahren!).
Die Jacobi-Matrix enthält die Werte der partiellen Ableitungen aller Zustandsgrößen nach allen Zustandsgrößen (=Wechselwirkung zwischen Zustandsgrößen).
Die Jacobi-Matrix wird erst neu berechnet, wenn dies auf Grund von Genauigkeitskriterien als erforderlich erkannt wird (z.B. mangelnde Konvergenz). Die Berechnung erfolgt in modernen Systemen symbolisch, in älteren noch über Abtastung des Modells.

Auswahl von Integrationsverfahren

In jedem kommerziellen Simulationsprogramm zur Dynamik-Simulation sind mehrere Integrationsverfahren implementiert, welche unterschiedliche Qualitäten besitzen. Der Nutzer sollte für seine konkrete Simulationsaufgabe ein Verfahren wählen, welches:

mit möglichst wenigen Stützstellen (Modellberechnungen) den zeitlichen Verlauf der Zustandsgrößen hinreichend genau ermittelt und
möglichst unempfindlich auf die Änderungen von Modellparametern reagiert (also stabil rechnet).

Die Anzahl der erforderlichen Modellberechnungen für ein konkretes Modell wird wesentlich von der maximalen Integrationsschrittweite dt bestimmt, mit welcher ein Verfahren noch stabil rechnet.
Um Verfahren in Hinblick auf die maximal zulässige Schrittweite dt vergleichen zu können, wird diese in Relation gesetzt zur kleinsten Zeitkonstante T bzw. zur Periodendauer T der höchsten Eigenfrequenz im Dynamikmodell (z.B. C·R, L/R, (L·C)^1/2, (m/c)^1/2 ).
Zusammen mit der Zahl der Stützstellen pro Integrationsschritt ergibt sich die Anzahl der erforderlichen Modellberechnungen n pro Zeitkonstante bzw. Periodendauer T.
Ein Kennzeichen für die Stabilität eines Verfahrens ist das Verhalten bei Überschreiten der maximal zulässigen Schrittweite. Günstig ist hierbei, wenn sich die Genauigkeit dabei nur stetig verschlechtert.

Explizite Verfahren
EULER                 : dtmax ≈ T/1000  n ≈ 1000/T (stark aufschaukelnd)
HEUN                  : dtmax ≈ T/10    n ≈   20/T (schwach aufschaukelnd)
Runge-Kutta 4.Ordn.   : dtmax ≈ T/5     n ≈   20/T (dämpfend, dann chaotisch)
Runge-Kutta-Fehlberg* : dtmax ≈ T/3     n ≈   18/T (dämpfend, dann chaotisch)
Implizite Verfahren (zusätzlicher Rechenaufwand für Erstellen der Jacobi-Matrix!)
Weighted Average      : dtmax ≈ T/200   n >  200/T (dämpfend, dann chaotisch)
Rosenbrock-Wanner*    : dtmax ≈ T/3     n >   18/T (stochastisch)
BDF-Verfahren*        : dtmax ≈ T/3     n >   18/T (stochastisch)

"Weighted Average" ist ein implizites Verfahren, welches mit einem Wichtungsfaktor in seiner Wirkung zwischen EULER- und HEUN-Verfahren eingestellt werden kann.
Die Verfahren mit geregelter Schrittweite (*) passen ihre aktuelle Schrittweite an die aktuell wirksame Zeitkonstante bzw. Periodendauer T an.
Auf den ersten Blick unterscheiden sich explizite und implizite Verfahren höherer Ordnung kaum im erforderlichen Berechnungsaufwand. Hier kommt jedoch ein wichtiger Unterschied zwischen beiden Verfahrensklassen zur Wirkung:
- Die Schrittweite dt eines expliziten Verfahrens richtet sich nur nach der kleinsten Zeitkonstante bzw. höchsten Eigenfrequenz, welche im Modell existiert. Dies gilt auch dann, wenn dieser kleine Wert für T aktuell nicht wirkt, weil z.B. ein harter mechanischer Anschlag aktuell nicht stattfindet! Das Modell muss also während der gesamten Simulationszeit unter Umständen mit extrem kleiner Schrittweite berechnet werden.
- Durch implizite Verfahren werden Zeitkonstanten und Eigenfrequenzen vernachlässigt, welche im aktuellen Modellverhalten nicht wirksam sind. Diese werden bei der Bildung der Jacobi-Matrix als "lineares Ersatzmodell" nicht berücksichtigt. Damit ist es möglich, steife DGL (mit extrem unterschiedlichen Eigenfrequenzen) überwiegend mit großen Integrationsschrittweiten zu berechnen.
Ein weiteres Problem sind die algebraischen Schleifen in den DGL, welche aus zusätzlichen Zwangsbedingungen zwischen Zustandsgrößen resultieren (z.B. Bewegungsgrößen in Getrieben):
- Aus obiger Liste sind dafür nur die BDF-Verfahren geeignet, welche auch als DASSL (Differential Algebraic System Solver) bekannt sind.
- Algebraische Schleifen sind auf Grund des Modellkonzepts in der Modellierungssprache Modelica die Regel. Deshalb wird bei Modelica-basierten Simulationssystemen meist ein BDF-Verfahren als Standardverfahren vorgegeben und benutzt.

Behandlung von Unstetigkeit als zeitdiskrete Ereignisse

Ereignisse im Simulationslauf

Spricht man von der Stetigkeit eines dynamischen Modells, so meint man damit nur die Stetigkeit der Zustandsgrößen-Funktionen:

Die Werte der Zustandsvariablen des Modells dürfen in stetigen Modellen zu keinem Zeitpunkt Sprünge machen.
Die Ableitungen der Zustandsgrößen dürfen in stetigen Modellen jedoch unstetig sein!

Die im Solver implementierten Verfahren der numerischen Integration liefern Lösungen nur innerhalb stetiger Bereiche dynamischer Modelle:

Die Grenze stetiger Modell-Bereiche ist durch sprunghafte Wertänderung von Zustandsgrößen gekennzeichnet.
Der Übergang in einen benachbarten stetigen Modell-Bereich muss durch eine sogenannte Ereignisbehandlung erfolgen, welche konsistente Anfangswerte für alle Zustandsgrößen berechnet, nachdem einzelne Zustandsgrößen sprunghaft ihren Wert geändert haben.
An der Grenze eines neuen stetigen Modell-Bereiches muss ein Neustart der numerischen Integration erfolgen (z.B. Erzeugen der "Historie" bei Mehrschrittverfahren).

Umgangssprachliche Definition von Ereignissen

In der Wikipedia kann man nachlesen, dass Ereignis vom neuhochdeutschen Wort eräugen abstammt:

Das Zustandekommen einer bestimmten Konstellation muss als Ereignis definiert (benannt) sein (meist vor dem Ereignis).
Es muss zusätzlich ein Beobachter existieren, damit Ereignisse bemerkt ("eräugt") werden.

Als Ereignis wird meist der Augenblick des Zustandekommens der besonderen Umstände gewertet. Es tritt dann zu einem diskreten Zeitpunkt schlagartig ein (keine Zwischenzustände!). Man kann zwei Formen von Ereignissen unterscheiden:

Prognostizierbare Ereignisse : Der Ereignis-Zeitpunkt ist durch Beobachtung der sich ändernden Umstände voraussagbar.
Zufällige Ereignisse : Wenn sich die Umstände aus Sicht des Beobachters stochastisch ändern, ist der Ereignis-Zeitpunkt erst nach dem Eintritt des Ereignisses bekannt.

Numerische Definition von Ereignissen

Die Unstetigkeiten eines Dynamik-Modells sind durch konkrete qualitative Merkmale gekennzeichnet (z.B. die Massen kollidieren). Damit solche Unstetigkeiten während der Simulation mit geringem numerischen Aufwand als Ereignisse erkannt werden können, wird die Definition eines Ereignisses auf einen einfachen Zahlenvergleich reduziert:

Erreichen eines Zustands während des Modellaufes, in dem eine sich zeitlich ändernde Größe den Wert einer Vergleichsgröße annimmt:
Als Beoabachter für jedes definierte Ereignis wird eine sogenannte Zero-Function genutzt. Diese registriert, bewertet und behandelt die Nulldurchgänge der Differenz von Ist- und Sollgröße.
Zero-Functions werden automatisch für alle unstetigen Funktionen (sign, abs, ...) und für die Testbedingung in Alternativ- und Schleifenanweisungen (if, for, ...) generiert.

Klassifizierung numerischer Ereignisse

Zufällige Ereignisse spielen nur eine Rolle bei der Fehlerbehandlung in Hinblick auf unzulässige Modell-Werte, wenn diese aus dem stochastischem Verhalten des Solvers oder eingespeisten Mess-Signalen resultieren:

Vermeidung mathematisch nicht definierter Operationen (z.B. Null-Division, negative Wurzel)
Vermeidung physikalisch-technisch ungültiger Modell-Werte (z.B. negative Masse, negative Kelvin-Absoluttemperatur), negativer Abstand)

Eine Ereignisbehandlung im Sinne der Überwindung von Modell-Unstetigkeiten wird hierfür im Normalfall nicht durchgeführt. Die Maßnahmen werden sich bei der Fehlerbehandlung auf das Setzen zulässiger Modell-Werte bzw. auf den Abbruch der Modellberechnung beschränken.

Bei den in den Dynamik-Modellen definierten Ereignissen handelt es sich überwiegend um prognostizierbare Ereignisse. Deren Ereignis-Zeitpunkt kann auf Grundlage der sich stetig ändernden Modellgrößen über eine begrenzte Zeitspanne hinweg vorausberechnet werden. Es wird dabei unterschieden zwischen zustands- und zeitabhängigen Ereignissen.

Zustandsabhängige Ereignisse:

werden durch Zerofunctions erkannt (ändern zum Ereigniszeitpunkt ihr Vorzeichen)
in unstetigen Funktionen (z.B. abs(x), sign(x))und bedingten Anweisungen (z.B. if..then..else)
Solver versucht, durch Prognose den Zeitpunkt möglichst genau zu treffen
Die kontinuierliche Integration wird kurz hinter dem tatsächlichen Ereigniszeitpunkt gestoppt (z.B. in SimulationX maximal dtDetect danach).
Zu diesem prognostizierten Zeitpunkt wird die Ereignisbehandlung gestartet

Zeitabhängige Ereignisse:

in zeitdiskreten Signalgliedern, z.B.
- AD-/DA-Wandler
- Abtast- & Halteglied
- diskrete Integratoren / Differentierer
- Z-Übertragungsfunktion
in Funktionen sample() bzw. delay() zum Setzen von Zeit-Ereignissen bzw. Signal-Verzögerung
Eintrag von zeitabhängigen Ereignissen in eine Ereigniswarteschlange - der Solver startet exakt zu diesen Zeitpunkten die Ereignisbehandlung

Ereignisbehandlung

Die Ereignisbehandlung soll einen neuen konsistenten Zustand aller Zustandsgrößen nach der Unstetigkeit herstellen:

Während der Ereignisbehandlung bleibt die Modellzeit konstant (time bzw. t).
Zu diesen Zeitpunkt wird eine so genannte Ereignis-Iteration durchgeführt.
Die Werte der elementaren Relationen (true bzw. false) bleiben während der kontinuierlichen Integration konstant, d.h. kurz vor einer Ereignis-Iteration wird noch mit den alten Werten gerechnet, obwohl sich die Relation schon geändert haben kann (Stopp nach dem wirklichen Ereigniszeitpunkt!)
Während der Ereignis-Iteration ändern die Relationen ihre Werte und danach werden für den gleichen Zeitpunkt neue konsistente Anfangswerte bestimmt. So erhält man zwei Werte eines Werteverlaufs an dem Zeitpunkt der Unstetigkeit: einen vor der Iteration und einen danach.
Mit einem noEvent-Operator kann man verhindern, dass ausgewählte Relationen Ereignisiterationen auslösen. Das kann man z.B. dazu benutzen, um den Definitionsbereich einer Funktion abzusichern:

y:= if noEvent(x>=0) then sqrt(x) else 0;

Der noEvent-Operator verhindert im Beispiel die automatische Generierung einer Zerofunction zur Ereignisbehandlung für die Testbedingung (x>=0). Ohne noEvent würde ein Fehler (Wurzel aus neg. Zahl) auftreten, weil kurz vor der Ereignis-Iteration der Wert von x schon negativ wird, obwohl x>=0 noch "true" bleibt.

Ein weiteres Beispiel aus der SimulationX-Hilfedatei zur die Vermeidung unerwünschter Ereignisbehandlung:

Als Sollwertverläufe für Positioniervorgänge werden gern Abschnitte von sin²-Funktionen verwendet:
Ein solches Signal kann man erzeugen mit dem Ausdruck:

if t<=0.5 then sin(pi*t)^2 else 1

Für die Relation wird automatisch ein Ereignis bei 0.5 s erzeugt, obwohl dort gar keine Unstetigkeit auftritt.
Der noEvent Operator führt hier zum Ausbleiben des Ereignisses und damit zu einer Beschleunigung der Rechnung, ohne dass damit Genauigkeitseinbußen verbunden sind:

if noEvent(t<=0.5) then sin(pi*t)^2 else 1

Beispiel: Extremwert-Sensor

Eine typische Modellierungsaufgabe ist die Erfassung des Extremwert-Betrages einer physikalischen Größe während des Simulationslaufes. Im Beispiel wird in einem SimulationX-Modell die Maximalgeschwindigkeit einer translatorischen Masse erfasst. Das Programmieren der Maximalwert-Erfassung erscheint auf den ersten Blick sehr einfach, wenn man zuvor den Extremwert vMax mit dem Anfangswert=0 versehen hat:

vMax := if abs(v) > vMax then abs(v) else vMax;

Leider gestaltet sich die Einbindung in ein Dynamik-Modell infolge des Wirkens der numerischen Integration und der Ereignisbehandlung ziemlich kompliziert:

Vor der Dynamik-Berechnung im Zeitbereich muss die Extremwert-Variable vMax einen sinnvollen Startwert erhalten (vMaxInit als Anfangswert von vMax).
Man muss gewährleisten, dass die "irreversible" Aktualisierung der aktuellen Extremwerte nur auf Basis "gültiger" Modellwerte erfolgt (last-Function):

Erläuterung:

Die Verwendung unstetiger Funktionen sollte nicht zu unnötigen Ereignisbehandlungen führen (noEvent-Operator).
Entnimmt man die zu überwachende Größe einem energetischen Konnektor, so muss man zur Vermeidung undefinierter Rückwirkungen die Flussgröße auf Null setzen (ctr1.F:=0):
Tritt der Extremwert nicht an einer bereits behandelten Unstetigkeitsstelle auf, so muss gezielt eine solche Testbedingung installiert werden! Nur so wird auch der Extremwert-Zeitpunkt berechnet (z.B. Nulldurchgang von a).
Infolge der Optimierung des Gleichungssystems vor der Simulation kann es passieren, dass in späteren Programm-Versionen solche zusätzlichen Test-Bedingungen "wegoptimiert" werden!
Sensor-Elemente von Modellbibliotheken liefern nicht nur "gültige" Signale. Das muss man ebenfalls durch Anwendung der last-Function berücksichtigen.

An diesem Beispiel sollte deutlich werden, dass die Berücksichtigung der Ereignisbehandlung bei der Modellbildung nicht einfach ist. Das Funktionieren von Tricks zur exakten Erfassung des Ereigniszeitpunkten ist stark abhängig vom "Optimierungspotential" des aktuell verwendeten Solvers.

Lineare und nichtlineare Systeme

Modellcharakter von Systemen

Der Begriff "System" wird inflationär mit verschiedenen Zusätzen benutzt, z.B.:

lineares bzw. nichtlineares System
chaotisches System
offenes, geschlossenes, abgeschlossenes System
emergentes System

Dies führt nicht unbedingt zu einem besseren Verständnis für das Wesen des Begriffs "System". Deshalb hier die Zusammenfassung der wesentlichen Merkmale, die man einem System zuordnet:

besitzt eine Abgrenzung zu seiner Umgebung;
enthält Komponenten, die Wirkung aufeinander auswirken;
es entsteht ein Systemverhalten, welches eine neue Qualität im Vergleich zu den Elementen aufweist (z.B. Spule und Kondensator ergibt einen Schwingkreis);
die resultierenden Systemgesetze engen den Freiheitsgrad der Komponenten ein! (Versklavung der Komponenten durch Abhängigkeiten).

Man muss unterscheiden zwischen

dem physikalischen System (real existierendes physikalisches Objekt) und
dem System-Konstrukt, welches als idealisiertes, vereinfachtes Abbild im Betrachter des physikalischen Systems entsteht.

Das System-Konstrukt ist somit ein Modell zum Verstehen eines als "System" betrachteten komplexen Gebildes.

Dynamische Systeme

Unter einem dynamischen System versteht man das mathematische Modell eines zeitabhängigen Prozesses, dessen Verlauf vom Anfangszustand abhängt:

Das System enthält Speicherelemente, d.h. es existiert ein Vektor von kontinuierlichen bzw. diskreten Zustandsgrößen Y(t).
Es gibt eine Überführungsfunktion T (=dynamisches Modell des Systems), welche die Änderung der Zustandsgrößen in einem Zeitintervall Δt beschreibt:

mit

Der Eingangsvektor u(τ) repräsentiert externe Signale, welche in das System eingespeist werden.

Man kann folgende Sonderformen von dynamischen Systemen unterscheiden:

deterministisches System: die Überführungsfunktion T ist eindeutig
stochastisches System: die Überführungsfunktion T ist mehrdeutig (Wahrscheinlichkeitsverteilung für das Erreichen der folgenden Zustände).
autonomes System: die Überführungsfunktion T ist unabhängig von der aktuellen Systemzeit t (d.h., das System entwickelt sich unabhängig vom Startzeitpunkt t₀).

Beispiele:

Ein Feder-Masse-Schwinger lässt sich als kontinuierliches, deterministisches und autonomes System beschreiben mit einer gewöhnlichen DGL beschrieben.
Biotope lassen sich als autonome, stochastische Systeme beschreiben (Übergangswahrscheinlichkeiten zwischen zeitlich aufeinanderfolgende Zustände).

Wann kann man ein System als lineares System behandeln

Für die Behandlung linearer Systeme existiert ein umfangreiches mathematisches Methoden-Inventar.
Das Verhalten linearer Systeme kann analytisch durch Anwendung der entsprechenden mathematischen Verfahren auf der Basis linearer Gleichungssysteme berechnet werden.

Beispiel "Theorie der Wechselstromschaltung":
Wenn die Schaltelemente konstante Parameter R, L, C besitzen und die Schaltung mit Sinusspannung gespeist wird, kann man durch Transformation der Differentialgleichung in den Bildbereich der komplexen Zahlen lineare Gleichungssysteme aufstellen und lösen.

Verallgemeinert gilt:
Sind Speicher- und Übertragungselemente konstant, so kann man Transformationen in einen Bildbereich durchführen, wenn sich der Eingangsvektor u(t) des Systems auf ein Spektrum harmonischer Schwingungen zurückführen lässt:

Sinus ⇒ Bildbereich der komplexen Zahlen,
Sprungfunktion und DIRAC-Impuls ⇒ Laplace-Transformation,
beliebige periodische Vorgänge ⇒ (Z-Transformation.

Wann muss man ein System als nichtlineares System behandeln

Wenn eines der folgenden drei Kriterien zutrifft, kann man ein System nicht mehr mit den Verfahren für lineare Systeme im Bildbereich berechnen:

Das System enthält speichernde Komponenten, deren Fassungsvermögen nicht konstant ist (z.B. Drehträgheit bei Eiskunstlauf-Pirouette)
oder
Das System enthält Übertragungselemente zwischen den Speichern, deren Übertragungsfunktion abhängig ist vom aktuellen Wert der übertragenen Größe (durch sich ändernde Werkstoffeigenschaften)
oder
Der Eingangsvektor u(t) ist aus Sicht des Systems zufällig.

Eine Transformation in einen Bildbereich und die analytische Lösung von linearen Gleichungssystemen ist dann nicht mehr möglich. In der Simulationspraxis hat dies im Normalfall nur geringe Bedeutung:

DGL beschreiben den Fluss von Energie, Stoff oder Information zwischen den Speicherkomponenten zeitkontinuierlicher dynamischer Systeme.
Auch für nichtlineare Systeme ist das Aufstellen der DGL meist einfach. Nur ihre analytische Lösung ist äußerst schwierig bis unmöglich.
"Zum Glück" können diese nichtlinearen dynamischen Modelle numerisch mit geeigneten Simulationssystemen genauso einfach berechnet werden. Allerdings sind dabei einige Besonderheiten zu beachten, welche im nächsten Abschnitt beschrieben werden.
Nur für Anwendungsfälle, in denen kurze Berechnungszeiten oberste Priorität besitzen, spielt die Möglichkeit der linearen Berechnung in einem Bildbereich zurzeit noch eine Rolle (z.B. bei Echtzeitfähigkeit von Modellen).

Numerische Stabilisierung unstetiger nichtlinearer Elemente

Die Verfahren der numerischen Integration funktionieren stabil auch für nichtlineare Systeme, wenn im zu berechnenden Modell alle erforderlichen zeitlichen Ableitungen der Zustandsgrößen einen stetigen Verlauf besitzen. Ein stetiger Verlauf (ohne Sprünge bzw. Polstellen) ist meist gewährleistet:

bei analytischer Formulierung der Zusammenhänge (z.B. physikalischer Effekte) ohne diskretisierende Operationen (wie Rundung, Vorzeichen, Fallunterscheidung) und ohne mögliche Divisionen durch Null (z.B. an Polstellen trigonometrischer Funktionen).
bei der Interpolation zwischen Stützstellen von Übertragungsfunktionen, welche nur in abgetasteter Form vorliegen (z.B. im Ergebnis von Messungen oder FEM-Simulationen), in Form von Polynomen oder Splines.

Unstetigkeiten in den Ableitungen führen meist zu Ereignisbehandlungen, welche den Simulationsfortschritt zumindest "ausbremsen" und bei Polstellen (Wert → ∞) meist zum abnormalen Abbruch der Simulation führen. Für einige in numerischer Hinsicht problematische Modell-Elemente werden nachfolgend Lösungen vorgeschlagen, wie Unstetigkeiten durch stetige Funktionen nachgebildet werden können, welche eine numerisch stabile Simulation ermöglichen.

Nichtlineare Energiespeicher

Das idealisierte, konzentrierte Element eines Energiespeichers wird beschrieben durch einen Parameterwert p (z.B. m, C, L):

Dieser Parameter p beschreibt vereinfacht, wie "träge" sich die Zustandsgröße ändert, wenn Energie zugeführt oder entzogen wird.
Für den Fall konstanter Speicherelemente (p=konstant) ergeben sich die bekannten einfachen Berechnungen der Zustandsgrößen:

Translatorische Bewegung einer Masse:
Rotatorische Bewegung einer Drehträgheit:
Strom durch Induktivität einer Spule:
Spannung einer elektrischen Kapazität:

Diese Zustandsgleichungen für konstante Energiespeicher p (d.h. mit dp/dt = 0) resultieren aus der sinngemäßen Anwendung des Impulserhaltungssatzes der klassischen Mechanik auf die unterschiedlichen physikalischen Partialsysteme:

Kraft verknüpft mit Impulsänderung:

Translatorische Bewegungen sind durch konstante Massen gekennzeichnet! Erst durch relativistische Effekte könnte es zu einer Masseänderung eines bewegten Körpers kommen. Masseänderungen durch Verdunstung oder Kondensation müssen wie Mehrkörpersysteme behandelt werden. Ein Kondensatteilchen verändert zwar die Masse des zu bewegenden Körpers, bewirkt aber keine Beschleunigung, da es seinen Impuls mitnimmt bzw. mitbringt.
Moment verknüpft mit Drehimpulsänderung:

Die rotatorische Drehträgheit ist eine Funktion der Massenverteilung in Bezug auf den Drehpunkt. Je mehr der konstanten Masse sich vom Drehpunkt entfernt, desto größer wird die Drehträgheit.
Induzierte Spannung verknüpft mit Induktionsflußänderung:

Die Induktivität der Spule eines Elektromagneten ist in erster Näherung eine Funktion des Luftspalts zwischen den beweglichen Teilen des Eisenkreises. Bei Berücksichtigung der Nichtlinearitäten des Eisenkreises ist die Induktivität auch noch abhängig von der aktuellen Flußdichte-Verteilung im Eisen, welche außer vom Luftspalt auch vom Spulenstrom bestimmt wird.
Strom in Kapazität verknüpft mit Ladungsänderung:

Die Kapazität eines kapazitiven Sensors ist eine Funktion des Elektrodenabstandes und der Größe der sich gegenüberstehenden Elektrodenflächen. Die Kapazität einer C-Diode ist eine Funktion der angelegten Spannung.

Falls es möglich ist, sollten die funktionellen Zusammenhänge zwischen dem Energiespeicher p und anderen physikalischen Größen in analytischer Form aus den physikalischen Effekten hergeleitet werden:

Gelingt dies nicht, kann man gemessene Kennlinienverläufe nutzen, welche z.B. als Kennfelder in Form diskreter Stützstellen (z.B. ASCII-Tabellen) bereitgestellt werden. Anstatt durch Messungen am realen Objekt kann man die erforderlichen Abhängigkeiten auch aus FEM-Rechnungen gewinnen (z.B. Abtastung eines 2D-Modell eines E-Magneten mit hinreichend vielen Stützstellen).
Es ist günstig, für solche Kennfelder möglichst jeweils eine geschlossene mathematische Ersatzfunktion zu identifizieren (z.B. p=f(x_i)).
Falls von dieser Funktion p=f(x_i) analytisch keine zeitliche Ableitung gebildet werden kann, kann man diese im Modell durch Anwendung der numerischen Differentation ermitteln. Es leidet darunter jedoch die numerische Stabilität der Modellberechnung, weil daraus Unstetigkeiten resultieren.
Besitzt man z.B. infolge variabler Energiespeicher eine weg-abhängige Funktion p(x) und benötigt dp(x)/dt, so vermeidet man die Probleme der numerischen Differentation mittels einen einfachen Tricks:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Formel Vermeidung num Differentation.gif

Die Bildung der partiellen Ableitung dp(x)/dx kann im Idealfall beim Vorliegen einer Funktion p(x) analytisch erfolgen.
Existiert p(x) nur als Kennfeld, so lässt sich daraus auf der Basis von Differenzenquotienten Δp(x)/Δx näherungsweise die partielle Ableitung dp(x)/dx ebenfalls als Kennfeld generieren, welches dann in der Simulation benutzt wird.
Die Geschwindigkeit v=dx/dt als Ursache der Wegänderung liegt im Dynamikmodell als stetige Zustandsgröße vor und wirkt damit stabilisierend auf die numerische Integration.

Viskose und Coulombsche Reibung

Bei der Betrachtung von Reibung kann man unterscheiden zwischen äußerer und innerer Reibung:

Äußere Reibung tritt als Reibkraft an den Berührungsflächen von Festkörpern auf und hemmt (Bewegungsreibung) bzw. verhindert (Ruhereibung) die Relativbewegung zwischen den Körpern.
Innere Reibung tritt innerhalb von Festkörpern, Flüssigkeiten und Gasen zwischen benachbarten Teilchen bei Verformungsvorgängen auf.

1. Innere (viskose) Reibung

Der Wert der inneren Reibungskräfte ist abhängig von der Relativgeschwindigkeit zwischen den betrachteten Komponenten.
Bei Festkörpern kann man mit Ausnahme von Gummimaterialien die innere Dämpfung meist vernachlässigen.
Hinweis: Der Fall der plastischen Verformung von Festkörpern ist stark nichtlinear und nicht durch Viskosität beschreibbar.
In speziell gestalteten Dämpfer-Elementen nutzt man die innere Reibung von Flüssigkeiten oder Gasen zur Erzeugung geschwindigkeitsproportionaler Relativkräfte. Daraus resultiert die nicht immer korrekte Bezeichnung "Viskose Reibung" für "Innere Reibung", da Dämpfer-Elemente die Hauptanwendung der inneren Reibung in der Technik darstellen.
Dämpfer-Modellelemente sind idealisierte Kraftelemente des mechanischen Partialsystems. Sie erzeugen Kraftkomponenten, welche der Relativbewegung zwischen Teilkörpern geschwindigkeitsabhängig entgegenwirken:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - viskose Daempfung.gif

für k > 0 gilt F1 = k * (v2-v1)

und es gilt immer: F2 = -F1

Dämpfungskräfte sind auch für variable "Dämpfungskonstanten" numerisch unkritisch, da es sich bei f(v) praktisch immer um eine stetige Funktion handelt.
Große Werte für die Dämpfungskonstante k entsprechen jedoch starken Wechselwirkungen zwischen den beteiligten Körpern. Dies erfordert z.B. bei Integrationsverfahren vom Runge-Kutta-Typ unzumutbar kleine Schrittweiten. Günstiger sind hier implizite Verfahren.

2. Äußere (coulombsche) Reibung

Nach der Art der Relativbewegung der Mechanismenteile zueinander unterscheidet man die folgenden Reibpaarungen:
Werden zwei Körper mit einer bestimmten Normalkraft F_N aneinandergepresst, so kann in dieser Reibpaarung eine Reibkraft F_R übertragen werden:
Das Verhältnis von Reibkraft zu Normalkraft wird Reibwert, Reibungszahl oder auch Reibungskoeffizient genannt:
Für die Reibpaarung wird ein idealisiertes konzentriertes Element eingeführt - das Reibelement:
Die Reibkraft F_R ist außer von der Normalkraft F_N auch von der Relativgeschwindigkeit Δv zwischen den Körpern abhängig. Sind beide Körper relativ zueinander in Ruhe (Δv=0), so spricht man von Ruhereibung, anderenfalls von Bewegungsreibung.
Ruhereibung: Dabei haften die Festkörperoberflächen mit ihren Rauhigkeiten aneinander (Haftreibung). Die durch die Reibpaarung miteinander wechselwirkenden Massen bleiben relativ zueinander in Ruhe, solange |F_Reib| < F_{Haft_max}.
- Die erzeugte Haftreibkraft wird als Reaktionskraft in Größe und Richtung nur durch die Summe der an den beiden Reibflächen angreifenden "inneren" Kräfte bestimmt (Actio=Reactio).
- Der Wert von µ_Haft läßt sich über den Grenzreibwinkel Ψ einer schiefen Ebene bestimmen:
Bewegungsreibung: Man unterscheidet drei Arten:

a) Festkörperreibung:

Hier bleiben die Oberflächen beider Körper infolge hoher Flächenkräfte dicht aneinander, so dass praktisch immer der Haftreibwert gilt ("Coulombsche Reibung" im ursprünglichen Sinne).

b) Flüssigkeitsreibung:

Sie tritt dann auf, wenn sich eine Gas- oder Flüssigkeitsschicht zwischen beiden Körpern befindet, z.B. bei einem Luftlager. Das Reibverhalten gleicht einer geschwindigkeitsproportionalen Dämpfung

c) Mischreibung:

Sie ist die am häufigsten auftretende Reibform. Ausgehend von dem Betrag des Haftreibungskoeffizienten nimmt der Reibwert µ_Reib mit zunehmendem Betrag der Geschwindigkeit ab und steigt anschließend wieder an.
Reibwert-Verläufe für Mischreibung zeigt das sogenannte Stribeck-Diagramm mit der Stribeck-Kurve.
Die exakte Bestimmung der Reibwert-Funktion ist nur schwer möglich! Reibwerte in Tabellenwerken sind recht grobe Näherungswerte.

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Reibwert-Tabellen.gif

Aus dem Reibwert-Verlauf µ_R=f(v) des Stribeck-Diagramms resultiert abhängig von der aktuell wirkenden Normalkraft F_N der folgende Kraftverlauf unter Berücksichtigung des wechselnden Vorzeichens der Relativgeschwindigkeit:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Reibkraft-Verlauf stetig.gif

Günstig für die Nachbildung des Reibwert-Verlaufs des Stribeck-Diagramms ist eine mathematische Funktion, welche durch die typischen Parameter dieses Reibwert-Verlaufs beschrieben wird, z.B.:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Reibwertverlauf qualitativ.gif

Man kann die Reibwert-Funktion aus 3 Bereichen zusammensetzen:

Haftreibung: Für Δv=0 wird my_Haft bereitgestellt.
Übergang von der Haft- zur Gleitreibung von Δv>0 bis zur Grenzgeschwindigkeit vGleitMin:
Dämpfungsanteil proportional Δv mit allmählicher Begrenzung bei vGleitMax:

Das Vorzeichen der Reibkraft muss in der Reibwert-Funktion nicht berücksichtigt werden. Es ergibt sich aus dem Vorzeichen der Relativ-Bewegung. Die vollständige Reibwert-Funktion für die Mischreibung lautet damit:

Besonderheiten der Haftreib-Modellierung:

Die Aufgabe besteht darin, eine solche Reibkraft zu ermitteln, dass zwei Massen im Zustand der Haftreibung numerisch weiterhin getrennt, aber hinreichend genau mit gleicher Beschleunigung nebenherlaufen:
Durch die Reibpaarung darf für den Fall Δv=0 nur soviel Kraft erzeugt werden, dass im Zusammenwirken mit den äußeren Kräften die reibenden Massen die gleiche Beschleunigung erfahren.
Der Kraftsprung bei Δv=0 widerspiegelt die Haftreibkraft, welche der aktuell angreifenden Kraft F entspricht und auf den Betrag der max. möglichen Haftreibkraft F_{Haft_max}=µ_Haft*F_N begrenzt ist, bei deren Überschreiten die Bewegung beginnt.
Der Vorzeichenwechsel der Reibkraft bei v=0 kann durch die Vorzeichenfunktion (-1, 0, +1) nachgebildet werden, wenn man die Start-/Stopp-Vorgänge beim Haft-/Gleitreib-Übergang mittels einer geeigneten Ereignisbehandlung organisiert. Diese Unstetigkeiten zusammen mit der Ereignisbehandlung führen jedoch zu numerischen Problemen!
Sprünge im Funktionsverlauf kann man grundsätzlich durch eine stetige Übergangsfunktion ersetzen. Da sich in mechanischen Systemen die meisten Unstetigkeiten auf Vorzeichenwechsel der Geschwindigkeit v oder des Abstandes x zurückführen lassen, spielt eine stetige Vorzeichenfunktion eine zentrale Rolle.
Der Übergang zwischen +1 und -1 wird hierbei durch einen Grenzwert für die verursachende Bewegungsgröße definiert (im Beispiel v). Bei einem hinreichend kleinem vGrenz unterscheidet sich die Wirkung der stetigen Signum-Funktion praktisch nicht mehr von der Original-Vorzeichenfunktion:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Signumfunktion stetig.gif

Wird der Haftreibungssprung im Nullpunkt durch eine steile Gerade zu ersetzen, die beim Erreichen der Grenzgeschwindigkeit vHaftMax kontinuierlich in den Wert my_Haft übergeht, so entspricht dies der Nachbildung einer zähen Dämpfungsschicht für den Haftreib-Fall:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Reibwertverlauf quantitativ.gif

Diese "extrem zähe Dämpfungsschicht" bewirkt eine Regelung der Reibkraft über die Abweichung der Geschwindigkeit von Null. Damit wird für den Kraftausgleich praktisch ein Regler mit Integralverhalten eingesetzt:

Die Anwendung der beschriebenen Prinzipe für ein translatorisches Reibelement zeigt der folgende Algorithmus im SimulationX TypeDesigner:

mit den zugehörigen Parameter und Variablen des Modell-Elements:

Grundlagen Simulation - Modellberechnung - nichtlin Elemente - Reibungselement Modellgroessen.gif

Hinweise:

Mittels der noEvent-Funktion wird konsequent die Ereignisbehandlung für unstetige Operationen deaktiviert.
Die Absolutwert-Bildung von x erfolgt teilweise als sqrt (x^2) zur Vermeidung der Ereignisbehandlung.
Zusätzlich zum Vorzeichenwechsel der Kraft bei v=0 wurde ein stetiger Übergang zwischen Festkörperreibung und Flüssigkeitsreibung bei v=vGleitMin in den Funktionsverlauf implementiert.

Wichtig:

Beim Design solcher stetiger Funktionsverläufe sollte man darauf achten, dass diese mittels der originalen Werkstoff-Parameter konfiguriert werden können (z.B. µ_Haft, µ_Gleit, k_Dämpfung)).
Die Konfiguration der Teilbereiche innerhalb des Funktionsverlaufes sollte durch augenscheinliche Positionen auf der Abszisse bestimmt werden (z.B. v_GleitMin, v_GleitMax).
Die "Weichheit" eines Funktionssprunges sollte ebenfalls durch einen anschaulichen Wert in Bezug auf die Absizzen-Skalierung definiert werden (z.B: v_HaftMax). Der Wert dieses Parameters muss zugeschnitten auf das aktuelle Modell als Kompromiss zwischen Berechnungsgeschwindigkeit und Genauigkeit des Modellverhaltens ermittelt werden (im Beispiel eine möglichst vernachlässigbare Restgeschwindigkeit im Haftreib-Zustand).

Mechanische Kontakt-Elemente

Kontakt-Elemente (auch als Anschlag-Elemente bezeichnet) werden in Modellen benötigt, um das Verhalten der Berührungsflächen bei Stoßvorgängen zwischen zwei Körpern nachzubilden:

Ein Stoß ist ein Vorgang, bei dem mindestens zwei Körper meist kurzzeitig Kraft aufeinander ausüben.
Im Verlauf eines Stoß-Vorganges ändern die Körper ihren Bewegungszustand, möglicherweise auch ihre Form und Zusammensetzung.
Für alle Stöße gelten der Impuls- und der Energie-Erhaltungssatz, wobei bei der Energie-Erhaltung auch die Umwandlung der Energieformen zu berücksichtigen ist.
Modelliert man ein Kontakt-Element als idealisierten starren Anschlag (d.h. mit der Zeitdauer=0 ohne Verformung der Berührungsflächen), so muss das Berührungsereignis erfasst werden und zu einem Initialisieren der resultierenden Körpergeschwindigkeiten entsprechend der Stoßgesetze führen (Siehe: "Stoß (Physik)" in Wikipedia)
Stoßvorgänge in mechanischen Anschlägen verringern häufig die Stabilität der numerischen Integration, wenn Sie die Ereignisbehandlung benötigen.

Im Folgenden werden elastische, stetige Ansätze für Kontakt-Elemente vorgestellt, welche unter Berücksichtigung der Materialeigenschaften den zeitlichen Verlauf von Stoßvorgängen einschließlich der wirkenden Kontaktkräfte im Modell nachbilden.

Unter einem "weichen" Stoßvorgang soll eine relativ weiche Berührung zweier Körperoberflächen (elastisch oder elastisch-plastisch) verstanden werden. Nach einem plastischen Stoß verbleiben beide Körper in Kontakt und bewegen sich gemeinsam weiter, bis sie durch äußere Kräfte getrennt werden.
Zur Behandlung von weichen Stoßvorgängen wird oft die Hertzsche Stoßtheorie herangezogen:

Die normierte Eindringtiefe Sn liegt bei Metall-Werkstoffen in der Größenordnung von 1 µm.

Da diese Stoßtheorie keine Energieverluste beinhaltet, bewirken die Stoßkräfte F periodische Schwingungen, was aber nicht dem realen Verhalten des Stoßes entspricht. Zur Ergänzung der Hertzschen Stoßtheorie muss deshalb noch eine Dämpfungskomponente eingeführt werden, die den Energieverlust des Stoßes berücksichtigt.
Die Federkraft F_F der elastischen Komponente ist von der Eindringtiefe Δs zwischen zwei stoßenden Körpern abhängig. Anstatt des Exponenten 0.5 sind auch andere Werte dafür möglich:
Für die Dämpfungskraft F_D unterscheidet man während der Berührung der Körperoberflächen zwei Fälle - die Kompressionsphase (Δv>0) und die Restitutionsphase (Δv<0):
Die Kontaktkraft F setzt sich dann aus den beiden Komponenten, Federkraft F_F und Dämpfungskraft F_D zusammen:

Die Parameter dieses Modellansatzes werden insbesondere bei metallischen Werkstoffen vor allem nach Kriterien der numerischen Stabilität und hinreichend geringer Eindringtiefe gewählt. Mit Exp=1/2 entspricht die Federkomponente F_F der Hertzschen Stoßtheorie.

Grundlage für die Einstellung der Dämpfungskonstante d im Kontakt-Element bildet das Lehrsche Dämpfungsmaß:

Bei der Masse m handelt es sich um die Gesamtmasse der am Stoß beteiligten Körper bzw. um die Masse des Körpers, welcher auf eine feste Oberfläche stößt.

D = 1 entspricht dem aperiodischem Grenzfall ohne Überschwingen.
D > 1 beschreibt ein stark gedämpftes Verhalten (Kriechfall).

Im Sinne der numerischen Stabilität ist der aperiodische Grenzfall günstig. Infolge der nichtlinearen Kraftkennlinien muss man in Abhängigkeit von Sn und exp einen günstigen Wert von D z.B. mittels einer Einflussanalyse ermitteln.
Diesen Modellansatz sollte nur für weiche Stoßvorgänge angewendet werden (d.h. für kleine Aufprallgeschwindigkeiten). Andernfalls kommt es beim Stoßvorgang zu einer unrealistisch großen Durchdringung s.

Modellierung eines elastischen Kontaktes in einem MKS-System (Autor: Sten Währisch):

Die im Rahmen von Forschungsprojekten am Institut für Feinwerktechnik und Elektronik-Design entwickelten Modellansätze zeigten, dass durch die Bestimmung der Steifigkeiten in allen drei Koordinatenrichtungen auch ein räumlicher elastischer Kontakt mit eindimensionalen Kraftelementen beschrieben werden kann, der ein Hertzsches Verhalten besitzt.
Die Steifigkeit des räumlichen Kontaktes ist ebenfalls abhängig von den Werkstoffparametern, der Geometrie und der Eindringtiefe.
Die Kontaktmodellierung mit Kraftelementen bleibt qualitativ natürlich hinter Oberflächen- oder Festkörpermodellen zurück. Dennoch liefert die Reduktion der Berührungsfläche auf einen Punktkontakt für spezielle Anwendungen (z. B. Hertzscher Kontakt mit Festkörperreibung) optimale Ergebnisse bei minimalem Aufwand.