Quelltext der Seite Grundlagen: Simulation Solver

<div align="center">[https://www.ifte.de/mitarbeiter/kamusella.html '''Autor: Dr.-Ing. Alfred Kamusella''']</div>

Simulationen mit numerischen Modellen beruhen im Wesentlichen auf der Berechnung von Modell-Ergebnissen mit Hilfe eines Gleichungslösers ([https://de.wikipedia.org/wiki/Solver 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 [https://de.wikipedia.org/wiki/SimulationX 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):<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_kenngroessen_simlauf.gif| ]]</div>
* 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:
[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_INIT-Abschnitt_-_beispiel_ode.gif|right]]

<span style="background-color:Yellow">  '''<u> INIT (&nbsp;''time''&nbsp;=&nbsp;''tStart''&nbsp;)</u>'''</span><br>
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: '''<u>O</u>'''rdinary '''<u>D</u>'''ifferential '''<u>E</u>'''quation). 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: [https://de.wikipedia.org/wiki/Differential-algebraische_Gleichung '''<u>D</u>'''ifferential '''<u>A</u>'''lgebraic '''<u>E</u>'''quation]). Hier können einige Zustandsgrößen nicht unabhängig voneinander gewählt werden:<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_INIT-Abschnitt_-_beispiel_dae.gif| ]]</div>
* 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.).

<span style="background-color:Yellow">'''<u>DYNAMIC (&nbsp;''tStart''&nbsp;≤&nbsp;''time''&nbsp;≤&nbsp;''tStop''&nbsp;)</u>'''</span><br>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.
 
<span style="background-color:Yellow">  '''<u> FINISH (&nbsp;''time''&nbsp;=&nbsp;''tStop''&nbsp;)</u>'''</span><br>
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 [https://de.wikipedia.org/wiki/Numerische_Integration numerischen Integration]:
* Zweck der Verfahren zur numerischen Integration ist die näherungsweise Berechnung  der Fläche ([https://de.wikipedia.org/wiki/Bestimmtes_Integral bestimmtes Integral]) unter einem Funktionsverlauf (Ableitung).
* Es existieren zwei grundsätzliche Möglichkeiten, um ein DGL-System zu notieren (hier am Beispiel des Pendels):
*# '''Ableitungsform''':<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_differentialform.gif| ]]</div>
*# '''Integralform''':<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_integralform.gif| ]]</div>
* 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 ===
[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_flaechenintegral.gif|right]]
'''Y(t)''' als Resultat einer Integral-Bildung kann geometrisch als Fläche unter dem Funktionsverlauf der zugehörigen Ableitung '''YP(t)''' interpretiert werden:<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_Wechselwirkung_Y_und_YP_-_integral.gif| ]]</div>
* '''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.: [https://de.wikipedia.org/wiki/Tangententrapezformel Trapez- /Tangentenformel] oder [https://de.wikipedia.org/wiki/Simpsonsche_Formel 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:<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_integrationskern.gif|Zentraler Integrationskern]] </div>

==== 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:<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_prinzip_num_integration.gif|.]] </div>
* Ausgehend von den Werten der Zustandsgrößen und Eingangssignale zum aktuellen Zeitpunkt '''t<Sub>n</Sub>''' 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<Sub>n</Sub>+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<Sub>n</Sub>''' können mit dem dynamischen Modell die aktuellen zeitlichen Ableitungen der Zustandsgrößen berechnet werden:  
<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_Ableitung_YP_berechnen.gif|.]] </div>
* 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. 


<span style="background-color:Yellow">'''<u> Explizite Verfahren </u>'''</span> 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<Sub>n</Sub>+dt''' als Anstieg der Tangente '''YP''' an die Zustandsgröße '''YP(t<Sub>n</Sub>)''' im Intervall '''dt''' dar:<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_explizite_verfahren.gif|.]] </div>
* Der prognostizierte Wert der Zustandsgrößen wird dann ebenfalls als hinreichend angenommen und ist Ausgangspunkt für die Berechnung des nächsten Zeitschrittes '''dt'''.  


<span style="background-color:Yellow">'''<u> Implizite Verfahren </u>'''</span> bewegen sich im Zeitbereich nicht so geradlinig vorwärts wie explizite Verfahren. Die Berechnung jedes Zeitschrittes '''dt''' ist durch eine Iterationsschleife gekennzeichnet:<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_implizite_verfahren.gif|.]] </div> 
* Unter Berücksichtigung des berechneten Anstiegs '''YP(t<Sub>n</Sub>)''' wird ein Anfangswert für die Zustandsgröße am Ende des Zeitschrittes '''dt''' '''Y(t<Sub>n</Sub>+dt)''' abgeschätzt (z.B. durch Extrapolation des bisherigen Signalverlauf '''Y(t)''').
* Der mit dem aktuellen Wert von '''Y(t<Sub>n</Sub>+dt)''' mittels des dynamischen Modells berechnete Anstieg '''YP(t<Sub>n</Sub>+dt)''' wird als Tangente an '''Y(t<Sub>n</Sub>+dt)''' angelegt und rückwärts zum aktuellen Zeitpunkt '''t<Sub>n</Sub>''' verfolgt.
* Ziel der anschließenden Iteration ist die Verbesserung des Schätzwertes für '''Y(t<Sub>n</Sub>+dt)''' solange, bis die rückwärts verfolgte Tangente '''YP''' den Wert von '''Y(t<Sub>n</Sub>)''' hinreichend genau trifft. 
* Ist dieses Ziel erreicht, werden die aktuellen Werte der Zustandsgrößen '''Y(t<Sub>n</Sub>+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.

<span style="background-color:Yellow">''' [https://de.wikipedia.org/wiki/Einschrittverfahren <u>Einschritt-Verfahren</u>] '''</span>
* "Betrachten" für die Prognose des nächsten Systemzustandes (Zeitpunkt='''t<sub>n</sub>+dt''') ausgehend von '''Y(t<sub>n</sub>)''' 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<sub>n</sub>+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).

<span style="background-color:Yellow">''' [https://de.wikipedia.org/wiki/Mehrschrittverfahren <u>Mehrschritt-Verfahren</u>] '''</span>
* 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 <span style="background-color:Yellow">''' [https://de.wikipedia.org/wiki/Explizites_Euler-Verfahren <u>EULER-Verfahren</u>] '''</span> ist das einfachste Verfahren:
* Man nimmt an, dass die Änderungsgeschwindigkeit '''YP''' der Zustandsgrößen '''Y''' über die Zeitspanne '''dt''' hinweg konstant bleibt:<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_Formel_Euler.gif|.]]</div>
* Infolge der Aufsummierung der Prognosefehler wird das EULER-Verfahren sehr schnell instabil!


Beim <span style="background-color:Yellow">''' [https://de.wikipedia.org/wiki/Heun-Verfahren <u>HEUN-Verfahren</u>] '''</span> wird ein Korrekturschritt angefügt:
* Mit der nach EULER prognostizierten Zustandsgröße Y<sub>1</sub>(t+dt) berechnet man die Ableitung YP(Y<sub>1</sub>(t+dt)).
* Dann nimmt man den Mittelwert beider Ableitungen als genauere Annahme für den Anstieg YP(t):<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_Formel_Heun.gif|.]]</div>
* 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 <span style="background-color:Yellow">''' [https://de.wikipedia.org/wiki/Runge-Kutta-Verfahren <u>Runge-Kutta-Verfahren</u>] '''</span>:
* Das EULER- und HEUN-Verfahren sind Runge-Kutta-Verfahren 1. bzw. 2. Ordnung.
* Je nach [https://de.wikipedia.org/wiki/Runge-Kutta-Verfahren#Beispiele Ordnung des konkreten Runge-Kutta-Verfahrens] wird der Funktionswert '''YP(t)''' auf der Basis mehrerer Stützstellen (Modellrechnungen) im betrachteten Zeitschritt "geschätzt":<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_Formel_allg_Runge-Kutta.gif|.]]</div>
* 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 [https://de.wikipedia.org/wiki/Runge-Kutta-Verfahren#Zeitschrittweitensteuerung:_Eingebettete_Verfahren 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 [https://de.wikipedia.org/wiki/BDF-Verfahren BDF-Verfahren] (Backward-Differential Formulas) handelt es sich um implizite Verfahren mit einstellbarer Ordnung und automatischer Schrittweitensteuerung:<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_BDF_als_implizites_Verfahren.gif|.]]</div>
* Auf Grundlage des bereits berechneten Verlaufs der Zustandsgröße '''Y(t)''' bis zum Zeitpunkt t erfolgt eine Startwert-Ermittlung '''Y<sub>o</sub>(t+dt)''' als Ausgangspunkt für die "exakte" Lösung '''Y(t+dt)''' je nach Ordnung des Verfahrens:<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_BDF_Startwert_Y0_fuer_naechsten_Schritt.gif|.]]</div>
* Danach erfolgt eine iterative Lösung des nichtlinearen Gleichungssystems zur Berechnung der Zustandsgrößen '''Y(t+dt)''' mit Abbruch bei der geforderten Genauigkeit.


Die <span style="background-color:Yellow">'''<u> BDF-Verfahren </u>'''</span> 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<sub>o</sub>(t+dt)''' extrapoliert (''Prädiktor''&nbsp;=&nbsp;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''&nbsp;=&nbsp;Verbesserer):<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_Korrektoriteration.gif|.]]</div>
* In jeder Korrektor-Iteration muss ein lineares Gleichungssystem gelöst werden:


<div align="center">'''res&nbsp;=&nbsp;J&nbsp;·&nbsp;(Y<sub>i-1</sub>&nbsp;-&nbsp;Y<sub>i</sub>)'''</div>
: mit  '''''res'''''&nbsp;=&nbsp;[https://de.wikipedia.org/wiki/Residuum_(Numerische_Mathematik) Residuen] / '''''J'''''&nbsp;=&nbsp;[https://de.wikipedia.org/wiki/Jacobimatrix Jacobimatrix] / '''''Y'''''&nbsp;=&nbsp;[https://de.wikipedia.org/wiki/Zustandsgröße Zustandsgrößen]
* Das Residuum widerspiegelt den momentanen Fehler in den Differentialgleichungen.
* Nach jeder Iteration wird die zu '''Y<sub>i</sub>''' gehörende Ableitung mittels Modellberechnung ermittelt und "rückwärts" überprüft, ob damit '''Y<sub>i-1</sub>''' 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)<sup>1/2</sup>, (m/c)<sup>1/2</sup> ).
* 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'''
 [https://de.wikipedia.org/wiki/Explizites_Euler-Verfahren EULER]                 : dtmax ≈ T/1000  n ≈ 1000/T (stark aufschaukelnd)
 [https://de.wikipedia.org/wiki/Heun-Verfahren HEUN]                  : dtmax ≈ T/10    n ≈   20/T (schwach aufschaukelnd)
 [https://de.wikipedia.org/wiki/Klassisches_Runge-Kutta-Verfahren Runge-Kutta 4.Ordn.]   : dtmax ≈ T/5     n ≈   20/T (dämpfend, dann chaotisch)
 [https://de.wikipedia.org/wiki/Runge-Kutta-Verfahren#Zeitschrittweitensteuerung:_Eingebettete_Verfahren Runge-Kutta-Fehlberg*] : dtmax ≈ T/3     n ≈   18/T (dämpfend, dann chaotisch)
 '''Implizite Verfahren''' ''(zusätzlicher Rechenaufwand für Erstellen der Jacobi-Matrix!)''
 [https://de.wikipedia.org/wiki/Implizites_Euler-Verfahren Weighted Average]      : dtmax ≈ T/200   n >  200/T (dämpfend, dann chaotisch)
 [https://de.wikipedia.org/wiki/Rosenbrock-Wanner-Verfahren Rosenbrock-Wanner*]    : dtmax ≈ T/3     n >   18/T (stochastisch)
 [https://de.wikipedia.org/wiki/BDF-Verfahren 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, [https://de.wikipedia.org/wiki/Steifes_Anfangswertproblem 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 [https://de.wikipedia.org/wiki/Differential-algebraische_Gleichung DASSL] (Differential Algebraic System Solver) bekannt sind.
** Algebraische Schleifen sind auf Grund des Modellkonzepts in der Modellierungssprache [https://de.wikipedia.org/wiki/Modelica 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 [https://de.wikipedia.org/wiki/Stetigkeit 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 [https://de.wikipedia.org/wiki/Ereignis 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:
# <span style="background-color:Yellow">'''<u> Prognostizierbare Ereignisse </u>'''</span>: Der Ereignis-Zeitpunkt ist durch Beobachtung der sich ändernden Umstände voraussagbar.
# <span style="background-color:Yellow">'''<u> Zufällige Ereignisse </u>'''</span>: 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:<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_num_ereignis.gif|.]]</div>
* 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 <u>nach</u> 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:  
<div align="center"> '''y:= if noEvent(x>=0) then sqrt(x) else 0;''' </div>
* 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:<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_num_ereignis_unerwuenscht.gif|.]] </div>
* Ein solches Signal kann man erzeugen mit dem Ausdruck:
<div align="center"> '''if t<=0.5 then sin(pi*t)^2 else 1''' </div>
* 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:
<div align="center"> '''if noEvent(t<=0.5) then sin(pi*t)^2 else 1''' </div>

==== 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:
<div align="center"> '''vMax := if abs(v) > vMax then abs(v) else vMax;''' </div>
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''):<div align="center">[[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_beispiel_extremwert_sensor.gif|.]]</div>
'''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 [https://de.wikipedia.org/wiki/System_(Begriffsklärung) "System"] wird inflationär mit verschiedenen Zusätzen benutzt, z.B.:
* [https://de.wikipedia.org/wiki/Lineares_System_(Systemtheorie) lineares] bzw. [https://de.wikipedia.org/wiki/Nichtlineares_System nichtlineares] System
* [https://de.wikipedia.org/wiki/Chaotisches_System chaotisches] System
* [https://de.wikipedia.org/wiki/Offenes_System offenes], [https://de.wikipedia.org/wiki/Geschlossenes_System geschlossenes], [https://de.wikipedia.org/wiki/Abgeschlossenes_System abgeschlossenes] System
* [https://de.wikipedia.org/wiki/Emergenz 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 [https://de.wikipedia.org/wiki/Physikalisches_System physikalischen System] (real existierendes physikalisches Objekt) und
* dem System-[https://de.wikipedia.org/wiki/Konstrukt 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:
<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_ueberfuehrungsfunktion.gif|.]] </div>
<div align="center"> mit </div>
<div align="center"> [[Datei:Grundlagen_Simulation_-_Modellberechnung_-_zeitkontinuierliche_Systeme_-_eingangsvektor.gif|.]] </div> 
* Der Eingangsvektor '''''u(τ)''''' repräsentiert externe Signale, welche in das System eingespeist werden.
Man kann folgende Sonderformen von dynamischen Systemen unterscheiden:
* [https://de.wikipedia.org/wiki/Systemeigenschaften#Determiniertheit deterministisches System]: die Überführungsfunktion T ist eindeutig
* [https://de.wikipedia.org/wiki/Stochastische_Systeme stochastisches System]: die Überführungsfunktion T ist mehrdeutig (Wahrscheinlichkeitsverteilung für das Erreichen der folgenden Zustände). 
* [https://de.wikipedia.org/wiki/Autonome_Differentialgleichung 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<sub>0</sub>).

'''Beispiele:'''
* Ein Feder-Masse-Schwinger lässt sich als kontinuierliches, deterministisches und autonomes System beschreiben mit einer gewöhnlichen DGL beschrieben.
* [https://de.wikipedia.org/wiki/Biotop 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 [https://de.wikipedia.org/wiki/Komplexe_Wechselstromrechnung "Theorie der Wechselstromschaltung"]:'''<br>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:'''<br>
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 [https://de.wikipedia.org/wiki/Harmonische_Schwingung#Harmonische_Schwingung harmonischer Schwingungen] zurückführen lässt:
* [https://de.wikipedia.org/wiki/Sinusschwingung '''Sinus'''] ⇒ [https://de.wikipedia.org/wiki/Komplexe_Wechselstromrechnung Bildbereich der komplexen Zahlen],
* [https://de.wikipedia.org/wiki/Sprungfunktion '''Sprungfunktion'''] und [https://de.wikipedia.org/wiki/Dirac-Impuls '''DIRAC-Impuls'''] ⇒ [https://de.wikipedia.org/wiki/Laplace-Transformation Laplace-Transformation],
* [https://de.wikipedia.org/wiki/Periodische_Funktion '''beliebige periodische Vorgänge'''] ⇒ ([https://de.wikipedia.org/wiki/Z-Transformation 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) <br>oder<br>
# Das System enthält Übertragungselemente zwischen den Speichern, deren Übertragungsfunktion abhängig ist vom aktuellen Wert der übertragenen Größe (durch sich ändernde Werkstoffeigenschaften)<br>oder<br>
# 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:''' [[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_konstante_Masse.gif|middle]]
# '''Rotatorische Bewegung einer Drehträgheit:''' [[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_konstante_Drehtraegheit.gif|middle]]
# '''Strom durch Induktivität einer Spule:''' [[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_konstante_Induktivitaet.gif|middle]]
# '''Spannung einer elektrischen Kapazität:''' [[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_konstante_Kapazitaet.gif|middle]]
Diese Zustandsgleichungen für konstante Energiespeicher '''p''' (d.h. mit '''dp/dt&nbsp;=&nbsp;0''') resultieren aus der sinngemäßen Anwendung des [https://de.wikipedia.org/wiki/Impuls ''Impulserhaltungssatzes der klassischen Mechanik''] auf die unterschiedlichen physikalischen Partialsysteme:
# '''<u>Kraft verknüpft mit Impulsänderung</u>:''' <br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Impulssatz_Masse.gif| ]] <br>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.<br>&nbsp;
# '''<u>Moment verknüpft mit Drehimpulsänderung</u>:''' <br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Impulssatz_Drehtraegheit.gif| ]] <br>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.<br>&nbsp;
# '''<u>Induzierte Spannung verknüpft mit Induktionsflußänderung</u>:''' <br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Impulssatz_Induktivitaet.gif| ]] <br>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.<br>&nbsp;
# '''<u>Strom in Kapazität verknüpft mit Ladungsänderung</u>:''' <br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Impulssatz_Kapazitaet.gif| ]] <br>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 [https://de.wikipedia.org/wiki/Kapazit%C3%A4tsdiode ''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<sub>i</sub>)).
* Falls von dieser Funktion '''p=f(x<sub>i</sub>)''' 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:<br>
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Vermeidung_num_Differentation.gif| ]]</div>
* 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. <u>Innere (viskose) Reibung'''</u>
* 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:
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_viskose_Daempfung.gif| ]]</div>
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_viskose_Daempfung.gif| ]]</div>
<div align="center">für '''k&nbsp;>&nbsp;0''' gilt '''F1&nbsp;=&nbsp;k&nbsp;*&nbsp;(v2-v1)'''</div>
<div align="center">und es gilt immer: '''F2&nbsp;=&nbsp;-F1'''</div>
* 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>Äußere (coulombsche) Reibung'''</u>
* Nach der Art der Relativbewegung der Mechanismenteile zueinander unterscheidet man die folgenden Reibpaarungen:<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Arten_von_Reibpaarungen.gif| ]]
* Werden zwei Körper mit einer bestimmten Normalkraft '''F<sub>N</sub>''' aneinandergepresst, so kann in dieser Reibpaarung eine Reibkraft '''F<sub>R</sub>''' übertragen werden:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Zusammenhang_zwischen_FN_und_FR.gif| ]]</div>
* Das Verhältnis von Reibkraft zu Normalkraft wird Reibwert, Reibungszahl oder auch Reibungskoeffizient genannt:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Reibwert.gif| ]]</div>
* Für die Reibpaarung wird ein idealisiertes konzentriertes Element eingeführt - das Reibelement:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibungselement.gif| ]]</div>
* Die Reibkraft '''F<sub>R</sub>''' ist außer von der Normalkraft '''F<sub>N</sub>''' 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<sub>Reib</sub>|&nbsp;<&nbsp;F<sub>Haft_max</sub>'''.
** 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 '''µ<sub>Haft</sub>''' läßt sich über den Grenzreibwinkel '''Ψ''' einer schiefen Ebene bestimmen:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibwert_mittels_Grenzwinkel_bestimmen.gif| ]]</div>
* '''Bewegungsreibung:''' Man unterscheidet drei Arten:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Arten_der_Bewegungsreibung.gif| ]]</div>
:a) <u>Festkörperreibung</u>:
:* 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) <u>Flüssigkeitsreibung</u>:
:* 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[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Mischreibung_im_Stribeck-Diagramm.gif|right]]
:c) <u>Mischreibung</u>:
:* Sie ist die am häufigsten auftretende Reibform. Ausgehend von dem Betrag des Haftreibungskoeffizienten nimmt der Reibwert '''µ<sub>Reib</sub>''' 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 [https://de.wikipedia.org/wiki/Stribeck-Kurve ''Stribeck-Kurve''].
:* Die exakte Bestimmung der Reibwert-Funktion ist nur schwer möglich! Reibwerte in Tabellenwerken sind recht grobe Näherungswerte.
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibwert-Tabellen.gif| ]]</div>
Aus dem Reibwert-Verlauf '''µ<sub>R</sub>=f(v)''' des Stribeck-Diagramms resultiert abhängig von der aktuell wirkenden Normalkraft '''F<sub>N</sub>''' der  folgende Kraftverlauf unter Berücksichtigung des wechselnden Vorzeichens der Relativgeschwindigkeit:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibkraft-Verlauf_stetig.gif| ]]</div>
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.:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibwertverlauf_qualitativ.gif| ]]</div>
* 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''':<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Uebergang_Haft-Gleit-Reibung.gif| ]]
# '''Dämpfungsanteil proportional Δv''' mit allmählicher Begrenzung bei '''vGleitMax''':<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Begrenzung_Reibdaempfung.gif| ]]
* 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:<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_vollstaendiger_Reibverlauf.gif| ]]
'''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:<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Haftreibung_zwischen_2_bewegten_Massen.gif| ]]</div>
* 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<sub>Haft_max</sub>=µ<sub>Haft</sub>*F<sub>N</sub>''' 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:
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Signumfunktion_stetig.gif| ]]</div>
* 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:
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibwertverlauf_quantitativ.gif| ]]</div>
* 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:
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Formel_Haftreibdaempfung_als_Integral-Regler.gif| ]]</div>
* Die Anwendung der beschriebenen Prinzipe für ein translatorisches Reibelement zeigt der folgende Algorithmus im ''SimulationX TypeDesigner'':<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibungselement_Algorithmus.gif| ]]</div>
: mit den zugehörigen Parameter und Variablen des Modell-Elements:
<div align="center">[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Reibungselement_Modellgroessen.gif| ]]</div>
'''''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. µ<sub>Haft</sub>, µ<sub>Gleit</sub>, k<sub>Dämpfung)</sub>).
* Die Konfiguration der Teilbereiche innerhalb des Funktionsverlaufes sollte durch augenscheinliche Positionen auf der Abszisse bestimmt werden (z.B. v<sub>GleitMin</sub>, v<sub>GleitMax</sub>).
* Die "Weichheit" eines Funktionssprunges sollte ebenfalls durch einen anschaulichen Wert in Bezug auf die Absizzen-Skalierung definiert werden (z.B: v<sub>HaftMax</sub>). 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: [https://de.wikipedia.org/wiki/Stoß_(Physik) '''"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 [https://de.wikipedia.org/wiki/Hertzsche_Pressung '''Hertzsche Stoßtheorie'''] herangezogen:<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Kontaktelement_Formel_nach_Hertzscher_Stosstheorie.gif| ]]
: Die normierte Eindringtiefe '''Sn''' liegt bei Metall-Werkstoffen in der Größenordnung von '''1&nbsp;µ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<sub>F</sub>''' 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:<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Kontaktelement_Formel_nach_Hertz_mit_freiem_Expon.gif| ]]
* Für die Dämpfungskraft '''F<sub>D</sub>''' unterscheidet man während der Berührung der Körperoberflächen zwei Fälle - die Kompressionsphase ('''Δv>0''') und die Restitutionsphase ('''Δv<0'''):<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Kontaktelement_Formel_daempfende_Komponente.gif| ]]
* Die Kontaktkraft '''F''' setzt sich dann aus den beiden Komponenten, Federkraft '''F<sub>F</sub>''' und Dämpfungskraft '''F<sub>D</sub>''' zusammen:<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Kontaktelement_Feder-_und_Daempfungskraft.gif| ]]
: 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<sub>F</sub>'''  der Hertzschen Stoßtheorie.
* Grundlage für die Einstellung der Dämpfungskonstante '''d''' im Kontakt-Element bildet das [https://de.wikipedia.org/wiki/Dämpfungsgrad '''Lehrsche Dämpfungsmaß]''':<br>[[Bild:Grundlagen_Simulation_-_Modellberechnung_-_nichtlin_Elemente_-_Kontaktelement_Formel_Lehrsches_Daempfungsmasz.gif| ]]
: 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&nbsp;=&nbsp;1''' entspricht dem aperiodischem Grenzfall ohne Überschwingen. 
:* '''D&nbsp;>&nbsp;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'''.
[https://www.optiyummy.de/images/Software_SimX_elastischer_Kontakt_im_MKS-System.pdf '''Modellierung eines elastischen Kontaktes in einem MKS-System''' (Autor: Sten Währisch)]:
* Die im Rahmen von Forschungsprojekten am [https://ifte.de ''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.