Thermodynamische Prozesse

Problemhintergrund

Als Beispiel in diesem Artikel dient eine typische Schnittstelle zu einem Hochdruck-Gasanschluss, wie sie bei verschiedenen Anwendungen wie der Steuerung von Aktoren (z. B. luftbetriebenen Kolben), bei pneumatischen Werkzeugen (z. B. pneumatischen Bohrern), in der Medizintechnik (Sauerstofftherapie, Beatmungsgeräte), industriellen Prozessanlagen usw. zum Einsatz kommt.

Das Modell (siehe Abbildung 1) besteht aus einem Druckregler und einem nachgeschalteten Steuerungsventil. Der Druckregler ist mit einer externen Druckquelle (tiefgestellte Ziffer 1) verbunden und sorgt für einen konstanten Eingangsdruck (tiefgestellte 2) an einem nachgeschalteten Steuerungsventil.

Der Druckregler hat an der Eingangsseite (Input) die Randbedingungen $p_1$ und $T_1$ und wir nehmen einen Druck von $p_3$ als Randbedingung für das Steuerungsventil an der Ausgangsseite (Output) an. Der Output des Druckreglers (tiefgestellte 2) ist nur idealisiert konstant; tatsächlich aber hängt er vom Durchfluss durch das System ab, daher ist er ein Zustand des Systems.

Abbildung 1: Thermodynamisches Modell mit einem Druckregler, einem Steuerungsventil und einem dazwischen befindlichen Volumen.

Zur Vereinfachung nehmen wir Folgendes an:

ideales Gas;
Druckregler und Steuerungsventil sind isenthalp;
es findet keine Gasmischung statt, d. h. es befindet sich nur ein einziges Gas im System.

Diese Vereinfachungen sind für viele Anwendungen sinnvoll. Wenn wir es jedoch mit grossen Druckunterschieden zu tun haben, kommt es bei einem realen Gas zu einer Temperaturänderung, wenn der Druck signifikant abnimmt (zum Beispiel aufgrund einer Öffnung). In diesem Fall wird durch die Berechnung mit einem idealen Gas diese Temperaturänderung nicht angezeigt. Wenn wir einen industriellen Prozess simulieren müssten, bei dem luftführende Rohre zur Konservierung mit Stickstoff befüllt werden müssen, läge eine Gasmischung vor, die folglich komplexere Modelle erfordern würde.

Die in den folgenden Abschnitten vorgestellten Modelle sind einfach, erfüllen aber den Zweck, ein thermodynamisches Netzwerk hinreichend genau zu simulieren, um bei der Entwicklung des Steuerungssystems von Nutzen zu sein.

Randbedingungen

Die Randbedingungen werden als konstant angenommen. An der Einlass-Seite haben wir $p_1$ und $T_1$ . An der Auslass-Seite brauchen wir lediglich den Druck $p_3$ , da kein Rückfluss gegeben ist ( $p_1$ > $p_3$ ).

Schritt 1: Mathematische Beschreibung des Modells

Zuerst betrachten wir jede Komponente einzeln und leiten ihre mathematische Darstellung her. Die Inputs, Outputs und Zustände werden bestimmt und die Gleichungen definiert.

Je nach Art der Komponente sind die Gleichungen entweder

gewöhnliche Differenzialgleichungen (GDGLs; engl.: ordinary differential equations, ODEs) erster Ordnung: für Komponenten mit stetigen Zuständen

$\frac{dy}{dt}=\ f(t,y,u)$

oder einfache explizite Gleichungen

$y\ =\ f(t,u)$

Darin ist t die Zeit, u der Input und y der Zustand bzw. Output.

Das Ventil

Bei gegebenem eingangs- und ausgangsseitigem Druck kommt es am Ventil zu einem Fluss (siehe Abbildung 2). Ein einfaches Ventil-Modell ohne modellierten Ventil-Aktor hat keine Zustände.

Abbildung 2: Pneumatisches Symbol eines Ventils mit den Inputs (Druck und Temperatur) des Modells und dem Massenfluss als Output.

Inputs:
- $p_2$ , $T_2$ , $p_3$ , $T_3$ , $Hub$
Outputs
- ${\dot{m}}_{Ventil}$
Zustände:
- keine

Die Massenfluss-Gleichung kann je nach Ventiltyp variieren, für Gase jedoch nimmt sie üblicherweise die folgende Form an:

${\dot{m}}_{Ventil}=\ 2\ast\ k\ast \sqrt{\frac{p_2 \ast \Delta p \ \ast \ p_3}{T_2 }} \ , wenn\ p_3>{\frac{p_2}{2}}$

${\dot{m}}_{Ventil}=\ k\ast\ p_2\ast\sqrt{\frac{\rho_2}{T_2}}\ ,\ sonst$

Mit = $k(Hub)$ ~ $K_V\left(Hub\right)$

Der Flusskoeffizient $K_V$ ist die standardisierte Art, um den Druckabfall mit dem Fluss über ein Ventil in Beziehung zu setzen. Er hängt von dem Ventilhub ab. Wenn $K_V$ nicht verfügbar ist oder die tatsächliche Kurve $K_V\left(Hub\right)$ von der der Konstruktion abweicht (z. B. aufgrund von Auswaschung), kann k aus den Messungen abgeleitet werden. Da das Ventil isenthalp und das Gas ideal ist, kommt es nicht zu einer Temperaturänderung, folglich ist $T_3$ = $T_2$ .

Der Druckregler
Der Druckregler (siehe Abbildung 3) reduziert einen (variierenden) Eingangsdruck $p_1$ auf einen Ausgangsdruck $p_2$ . Idealerweise wird der Ausgangsdruck auf den Sollwert-Druck geregelt. Der Druckregler besteht aus einer Feder, die das Ventil öffnet und der geregelte Druck erzeugt eine Gegenkraft, die das Ventil bei Drücken oberhalb des Sollwerts schliesst.

Abbildung 3: Pneumatisches Symbol eines Druckreglers mit den Inputs (Druck und Temperatur) des Modells und dem Massenfluss als Output.

Inputs:
- $p_1$ , $T_1$ , $p_2$ , $T_2$
Outputs:
- ${\dot{m}}_{Ventil}$
Zustände:
- keine

Ein einfacher Ansatz zur Modellierung des Druckreglers besteht in der Anwendung derselben Gleichung wie für das Ventil, allerdings mit $k\ =\ k\left(p2\right)$ (siehe Abbildung 4):

${\dot{m}}_{Reg}=\ 2\ast k(p_2)\ast \ \sqrt{\frac{p_1 * \Delta p * p_2}{T_1}}} \ , \ wenn \ p_2 > \frac{p_1}{2}$

${\dot{m}}_{Reg}=\ k\left(p_2\right)\ast p_1\ast\sqrt{\frac{\rho_1}{T_1}}\ ,\ sonst$

Abbildung 4: Druckabhängiger Flusskoeffizient für das Ventil-Modell.

Für einen Ausgangsdruck über dem Druckregler-Sollwert ist der Fluss gleich 0. Bei einem Druck unter einem Schwellenwert ist der Druckregler vollständig geöffnet. Dazwischen wird der Flusskoeffizient interpoliert. Dieser geringere Schwellenwert sollte nicht zu nahe am Sollwert-Druck sein, andernfalls kann es schwierig sein für einen GDGL-Löser, eine stabile und konvergierende Lösung zu berechnen.

Das Volumen

Zwischen dem Ventil und dem Druckregler befindet sich eine Rohrleitung gewisser Länge mit einem Volumen, die modelliert werden muss. Es nimmt einen Fluss ${\dot{m}}_{Reg}$ aus dem Druckregler auf und hat einen Ausfluss ${\dot{m}}_{Ventil}$ durch das Steuerungsventil (unter der Annahme, dass das Ventil offen ist und der Druck hinter dem Ventil geringer ist). Zusätzlich kann eine Energieübertragung $\dot{Q}$ in die bzw. aus der Umgebung eingeschlossen werden.

Abbildung 5: Volumen mit den Modell-Inputs Massenfluss, Temperaturen und Wärmestrom sowie den Zuständen Druck und Temperatur.

Inputs:
- ${\dot{m}}_{Reg}$ , $T_{Reg}$ , ${\dot{m}}_{Ventil}$ , $\dot{Q}$
Outputs:
- $p_2$ , $T_2$
Zustände:
- $p_2$ , $T_2$

Statt der Temperatur kann die Enthalpie als ein Zustand eingesetzt werden (zum Beispiel wenn das Gas ein Dampf ist und daher Dampftabellen verwendet werden müssen). Mit der Gleichung eines idealen Gases, der Massenerhaltung und der Energieerhaltung können wir die GDGLs herleiten:

$\dot{T}}_2=\frac{1}{m\ast c_V}\left(\left[-{\dot{m}}_{Reg}\ast c_V-{\dot{m}}_{Ventil}\ast R\right]\ast T_2+{\dot{m}}_{Reg}\ast T_{Reg}\ast c_p-\dot{Q}\right)$

${\dot{p}}_2=\ \frac{R}{V}\left(\left[{\dot{m}}_{in}-{\dot{m}}_{out}\right]\ast T_2+\ m\ast{\dot{T}}_2\right)$

Die Wärmeübertragung $\dot{Q}$ kann, falls erforderlich, zum Beispiel mit der natürlichen Konvektion, der forcierten Konvektion oder auch Strahlung (bei sehr hohen Temperaturen) modelliert werden.

Schritt 2: Das System lösbar machen

Die Gleichung für das Ventil (und auch für den Druckregler) zur Berechnung des Massenflusses

${\dot{m}}_{Ventil}=\ 2\ast\ k\ast \sqrt{\frac{p_2 \ast \Delta p \ \ast \ p_3}{T_2 }}}$

ergibt ein Problem für die Annäherung von $p_3$ an $p_2$ . Die Ableitung geht gegen unendlich

$\lim\underset{p_3\rightarrow p_2} \ {\left(\frac{d{\dot{m}}_{Ventil}}{dp_3}\right)\rightarrow\pm\infty}$

was bedeutet, dass eine kleine Änderung des $\Delta p$ im Bereich von 0 zu einer proportional grossen Flussänderung führt. Bei der Berechnung der Jacobi-Matrix für die erweiterten GDGL-Löser kann dies zu Stabilitätsproblemen führen.

Um dieses Problem zu umgehen, kann die Gleichung für kleines $\Delta p$ linearisiert werden:

${\dot{m}}_{Ventil}=\ 2\ast k\ast\sqrt{\frac{p_2 \ast \Delta p \ast p_3}{T_2}} , \ für \ \Delta p \ge \Delta p_m_i_n$

${\dot{m}}_{Ventil}=\ 2\ast k\ast\sqrt{\frac{p_2\ast p_3}{T_2}} \ast c \ast \Delta p , \ für \ \Delta p < \Delta p_m_i_n$

mit $c$ so, dass ein kontinuierlicher Übergang zwischen den beiden Gleichungen besteht.
Andere Ansätze mit Polynomen höherer Ordnung sind ebenfalls möglich. So ermöglicht zum Beispiel ein Polynom zweiter Ordnung nicht nur einen kontinuierlichen Übergang, sondern auch einen nahtlosen Übergang bei $\Delta p$ = $\Delta p_m_i_n$ .

Schritt 3: Die Modelle vereinen

Damit haben wir nun einen Satz von Gleichungen, die gelöst werden können. Das Ventil und der Druckregler liefern Massenflüsse, basierend auf den Zuständen und den Randbedingungen für das Volumen; aus dem Volumen berechnen sich die Ableitungen der Zustandsgrössen und der GDGL-Löser kann dann die neuen Zustände auf der Basis der abgeleiteten Werte berechnen. Abbildung 6 zeigt den Arbeitsablauf zur Lösung des Systems.

Abbildung 6: Arbeitsablauf zur Lösung des Systems.

Schritte:

Der Löser berechnet die aktuellen Zustände $T_2$ und $p_2$ .
Die Massenflüsse ${\dot{m}}_{Reg}\left(p_1,T_1,p_2\right)$ und ${\dot{m}}_{Ventil}\left(p_2,T_2,p_3,\ Hub\right)$ werden berechnet.
Die Ableitungen $T_2$ und $p_2$ werden berechnet.

Dies kann zum Beispiel in Simulink [1] umgesetzt werden, siehe Abbildung 7. Darin ist das System mit einem Euler-Vorwärtslöser dargestellt.

Abbildung 7: Simulink-Umsetzung des thermodynamischen Systems mit dem Euler-Vorwärtslöser.

Schritt 4: Das System lösen

Es gibt zwei Arten von GDGL-Lösern:

Löser mit fester Schrittweite
Löser mit variabler Schrittweite

Bei den Lösern mit fester Schrittweite wird eine konstante Schrittweite angenommen und der nächste Schritt basierend auf dem aktuellen und (abhängig von der Methode) vorherigen oder Zwischenschritten kalkuliert. Bei den Lösern mit variabler Schrittweite auf der anderen Seite wird die Schrittweite angepasst, um sicherzustellen, dass der Fehler pro Schritt nicht die relativen und absoluten Grenzwerte überschreitet. Diese Löser sind stabiler, da die Schrittweite unter instationären Bedingungen reduziert ist; das bedeutet aber auch, dass die Berechnungszeit in solchen Situationen erhöht ist. Die erforderliche Zeitspanne für die Berechnungen kann daher in Abhängigkeit von dem aktuellen Systemzustand stark variieren.

Für einfache Anwendungen kann ein Löser mit fester Schrittweite ausreichend sein. Für anspruchsvollere System jedoch ist ein Löser mit variabler Schrittweite häufig die bessere Wahl. In den nächsten Abschnitten werden einige verschiedene Löser untersucht, um ihre Leistung und Stabilität bei gegebenen verschiedenen Systemparametern zu zeigen.

Löser mit fester Schrittweite

Im Folgenden vergleichen wir drei Methoden:

Euler-Vorwärtsmethode (einfachste Runge-Kutta-Methode)
Heun-Verfahren (zweistufige Runge-Kutta-Methode)
dreistufige Runge-Kutta-Methode

Für die Anfangsbedingungen nahmen wir $p_{2,0}=\ 1\ bar$ und $T_{2,0} =\ 25\ ^{\circ}C$ an und als Randbedingung haben wir $p_1=\ 5\ bar$ , $T_1=\ 25\ ^{\circ}C$ und $p_3=\ 1\ bar$ . Für das Ventil nehmen wir einen konstanten Hub an.

Wir werden folgende Beispiele untersuchen:

Beispielnummer	Volumen	Schrittzeit
1	3 cm³	5 ms
2	2 cm³	5 ms
3	1.6 cm³	5 ms
4	1.6 cm³	1 ms

Beispiel 1
Beim ersten Beispiel setzen wir $V\ =\ 3 \ {cm}^3$ und $dt\ =\ 5\ ms$ . Das Ergebnis ist in Abbildung 8 dargestellt. Das Heun-Verfahren und die 3-stufige Runge-Kutta-Methode führen zu einem nahezu identischen Ergebnis sowohl für den Druck als auch die Temperatur, wohingegen die Temperatur bei der Euler-Vorwärtsmethode während des instationären Verhaltens ein wenig höher ist. Die Euler-Vorwärtsmethode ist die einfachste, aber auch die genaueste dieser drei Methoden.

Abbildung 8: Gelöstes System, mit Lösern mit fester Schrittweite, V = 3 cm³ und dt = 5 ms

Beispiel 2
Beim zweiten Beispiel setzen wir $V\ =\ 2 \ {cm}^3$ und $dt\ =\ 5\ ms$ . Das Ergebnis ist in Abbildung 9 dargestellt. Mit dem kleineren Volumen kommt es bei der Euler-Vorwärtsmethode anfangs zu einer Oszillation. Alle drei Methoden sind stabil und konvergieren.

Abbildung 9: Gelöstes System, mit Lösern mit fester Schrittweite. V = 2 cm³ und dt = 5 ms

Beispiel 3

Beim dritten Beispiel setzen wir $V\ =\ 1,6 \ {cm}^3$ und $dt\ =\ 5\ ms$ . Das Ergebnis ist in Abbildung 10 dargestellt. Bei dem noch kleineren Volumen verhalten sich Euler-Vorwärts- und 3-stufige Runge-Kutta-Methode weiterhin stabil, allerdings konvergieren sie nicht mehr. Das Heun-Verfahren scheint gute Ergebnisse zu liefern, da die Lösung stabil und konvergierend ist; bei Betrachtung der Temperatur bei 0,1 Sekunden ist jedoch offensichtlich, dass das Ergebnis nicht korrekt sein kann, da die Temperatur gegen die Eingangs-Randbedingungstemperatur von 25 °C konvergieren sollte.

Abbildung 10: Gelöstes System, mit Lösern mit fester Schrittweite. V = 1,6 cm³ und dt = 5 ms

Beispiel 4

Beim vierten und letzten Beispiel setzen wir $V\ =\ 1,6 \ {cm}^3$ und $dt\ =\ 1\ ms$ . Das Ergebnis ist in Abbildung 11 dargestellt. Mit der reduzierten Schrittweite konvergiert die Lösung nun wieder gegen eine stabile Lösung.

Abbildung 11: Gelöstes System, mit Lösern mit fester Schrittweite. V = 1,6 cm³ und dt = 1 ms

Die Anzahl der Calls zur Berechnung der Ableitung hängt von der Methode ab, sie wächst aber bei allen Methoden mit der Zeitdauer, die simuliert wird, linear an. Dies ist ein Nachteil dieser einfachen Methoden. Auch wenn die Systemzustände sich nicht mehr viel ändern, wird die Schrittweite nicht erhöht.

	Euler-Vorwärts	Heun	dreistufige Runge-Kutta
0.1 s	100	200	300
1 s	1’000	2’000	3’000
10 s	10’000	20’000	30’000
100 s	100’000	200’000	300’000

Löser mit variabler Schrittweite

Wir vergleichen zwei verschiedene Löser

ode45
lsode

Wir nehmen dieselben Randbedingungen wie bei den Lösern mit fester Schrittweite an. Das Volumen wird auf $V\ =\ 1,6 \ {cm}^3$ gesetzt.

ode45

Der ode45-Löser ist ein Standard-GDGL-Löser einer linearen Mehrschritt-Methode in Matlab [2]. Im Gegensatz zu den im vorigen Abschnitt vorgestellten Lösern mit fester Schrittweite speichert er die Informationen aus den vorherigen Schritten, um mehr Effizienz bei der Lösung des GDGL-Systems zu erzielen.

Das Simulationsergebnis des thermodynamischen Systems ist in Abbildung 12 dargestellt. Am Anfang beträgt die Schrittweite ungefähr 1 ms, sie wird anschliessend aber bis auf etwa 5 ms erhöht.

Abbildung 12: Mit dem ode45-Löser in Simulink gelöstes System. V = 1,6 cm³.

Die folgende Tabelle listet die Anzahl der Calls zur Berechnung der Ableitung auf.

	ode45
0.1 sec	289
1 sec	1’465
10 sec	13’111
100 sec	129’589

lsode – Livermore-Löser für GDGLs

Der Livermore-Löser für GDGLs [3] ist ein Löser für steife (Rückwärtsdifferenzierungs-Formel; engl. Backward Differentiation Formula, BDF) und unsteife (Adams-Methoden) Probleme mit einer Jacobi-Matrix, die vom Anwender eingegeben oder intern erzeugt werden kann. Der Löser liefert Lösungen an exakt definierten Zeitpunkten und speichert Informationen aus vorherigen Schritten für die nächsten Schritte.

Bei der Lösung eines GDGL-Systems als Teil einer grösseren Simulationsumgebung ist es üblicherweise erforderlich, das System schrittweise zu lösen und Lösungen für einen gegebenen Zeitpunkt zu haben. In einer HIL-/SIL-Umgebung gibt es einen diskreten Regler, der die Eingangsgrössen (Inputs) zum thermodynamischen System steuert (z. B. Ventile). Dieser Regler hat eine feste Schrittweite und alle x ms liest er die Inputs und stellt die Outputs ein. Daher muss der thermodynamische Prozess auch alle x ms die Ergebnisse ausgeben. Der lsode-Algorithmus kann diese Funktionalität bereitstellen. Für diese Simulation wurde das Ventil bei einem konstanten Hub gehalten, sodass das Ergebnis mit den anderen Simulationen verglichen werden kann.

Die Simulation des thermodynamischen Systems mit dem lsode-Löser kann in Abbildung 13 betrachtet werden. Die Simulation wurde in MS Visual Studio 2015 programmiert und gelöst. Die Simulation wurde mit einer Schrittweite von 5 ms durchgeführt. Der Löser könnte allerdings auch mehrere Teilschritte ausführen, die nicht gezeigt werden. Am Anfang sind viele Teilschritte erforderlich, um ein Ergebnis zu erhalten, das die relativen und absoluten Toleranzen einhält; gegen Ende jedoch kann die nächste Iteration in nur einem Schritt berechnet werden. Das Ergebnis mag etwas grob aussehen, weil die Teilschritte in der Abbildung nicht dargestellt sind.

Abbildung 13: Mit dem lsode-Löser gelöstes System. V = 1,6 cm³, Schrittweite = 5 ms.

Die folgende Tabelle listet die Anzahl der Calls zur Berechnung der Ableitung auf.

	lsode
0.1 s	221
1 s	281
10 s	302
100 s	306

Wenn die Schrittweite auf 1 ms gesenkt wird (siehe Abbildung 14), ähnelt das Ergebnis mehr dem Ergebnis des Systems, das mit dem ode45-Löser wie in Abbildung 12 gezeigt gelöst wurde, wo am Anfang viele Schritte mit einer kleinen Schrittweite gemacht wurden.

Abbildung 14: Mit dem lsode-Löser gelöstes System. V = 1,6 cm³, Schrittweite = 1 ms.

Vergleich

Ein Vergleich der verschiedenen Lösungen ist in Abbildung 15 dargestellt. Die Lösung des lsode-Lösers erscheint im Vergleich zu den anderen Lösungen sehr grob, dies liegt jedoch einfach daran, dass die Lösung nur für die erforderliche Schrittweite von 5 ms abgebildet ist, d. h. es wird nur die alle 5 ms erhaltene Lösung dargestellt, nicht jedoch die Lösung bei allen Teilschritten.

Abbildung 15: Vergleich der mit verschiedenen Lösern erhaltenen Lösungen des Systems.

Die folgende Tabelle listet die Anzahl der Calls zur Berechnung der Ableitung für jeden Löser auf.

	Euler-Vorwärts	Heun	3-stufige RK	ode45	lsode
0.1 s	100	200	300	289	221
1 s	1’000	2’000	3’000	1’465	281
10 s	10’000	20’000	30’000	13’111	302
100 s	100’000	200’000	300’000	129’589	306

[1]	„Simulink 9.2,“ The Mathworks Inc., Natick, Massachusetts, 2018.
[2]	„Matlab 9.5 (R2018b),“ The MathWorks Inc., Natick, Massachusetts, 2018.
[3]	A. Hindmarsh, „Ordinary Differential Equation Solvers,“ 13 02 2008. [Online]. Available: https://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html. [Accessed 04 08 2021].