Modellierung thermodynamischer Prozesse in Gassystemen

Eine Übersicht über die numerischen Gleichungslöser für Modelle thermodynamischer Prozesse auf der Grundlage eines Hochdruck-Gassystems

Die frühzeitige Systemverifizierung ist der Schlüssel, um unvorhergesehene Probleme während des Entwicklungsprozesses zu vermeiden. Aus diesem Grund wendet IMT mathematische Methoden nicht nur an, um die Umsetzbarkeit komplexer Systeme zu prüfen, sondern auch um mögliche Konstruktionsprobleme in einer frühen Phase zu identifizieren. Das Verfahren umfasst üblicherweise, das System in Prozessmodellen mathematisch zu beschreiben und anschliessend anhand von Simulationen zu überprüfen.

Eine mathematische Darstellung des Systems zu haben, ermöglicht es IMT auch, den Regler-Entwicklungsprozess zu beschleunigen. Der Regler wird dazu emuliert und mit dem Prozessmodell verbunden. Das bietet die Möglichkeit, an dem Regler zu arbeiten, ohne ihn für das Gerät anfertigen zu müssen; darüber hinaus können so mit der Auflösung von nur einem Schritt Erkenntnisse gewonnen werden, wie der Regler sich verhält.

Im ersten Teil des Beitrags betrachten wir die Modellierung einiger typischer Komponenten für Hochdruck-Gasgeräte und leiten die Gleichungen zur Beschreibung der Systemdynamik her. Im zweiten Teil werden wird dann verschiedene Gleichungslöser vorstellen und Beispiele für unterschiedliche, von verschiedenen Parametern abhängende Lösungen geben, um ihre Vor- und Nachteile aufzuzeigen.

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.

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.

• 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.

• 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 * 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 \)

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.

• 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 ≥ \Delta p_{min} \)

\( {\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_{min} \)

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_{min} \).

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.

Schritte:

1. Der Löser berechnet die aktuellen Zustände \( T_2 \) und \( p_2 \).
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.
3. Die Ableitungen \( T_2 \) und \( p_2 \) werden berechnet.

Dies kann zum Beispiel in Simulink („Simulink 9.2,“ The Mathworks Inc., Natick, Massachusetts, 2018) umgesetzt werden, siehe Abbildung 7. Darin ist das System mit einem Euler-Vorwärtslöser dargestellt.

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:

1. Euler-Vorwärtsmethode (einfachste Runge-Kutta-Methode)
2. Heun-Verfahren (zweistufige Runge-Kutta-Methode)
3. 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:

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.

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.

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.

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.

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.

Löser mit variabler Schrittweite

Wir vergleichen zwei verschiedene Löser

1. ode45
2. 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 („MATLAB 9.5 (R2018b),“ The MathWorks Inc., Natick, Massachusetts, 2018.). 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.

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

lsode – Livermore-Löser für GDGLs

Der Livermore-Löser für GDGLs (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].) 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.

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

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.

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.

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

In diesem Artikel haben wir zunächst die Schritte der Modellierung eines einfachen thermodynamischen Systems betrachtet. Die Modellierung führt zu mehreren Gleichungen, die auf einfache Weise die Berechnung der Zustandsableitungen ermöglichen.

In einem zweiten Schritt haben wir mehrere Gleichungslöser betrachtet – Löser mit fester Schrittweite und Löser mit variabler Schrittweite. Der Vergleich zeigt eindeutig den Vorteil der Löser mit variabler Schrittweite, da sie die Schrittweite an das Systemverhalten anpassen können. Darüber hinaus zeigt lsode, der Löser mit variabler Schrittweite, einen klaren Vorteil gegenüber ode45, ebenfalls einem Löser mit variabler Schrittweite, da er signifikant weniger Schritte erfordert und somit das System wesentlich schneller löst. Der ode45-Löser ist ausserdem an MATLAB gebunden, jedoch können die Löser mit fester Schrittweite in einer beliebigen Sprache programmiert und der lsode-Löser kann als Bibliothek in die meisten geläufigen Sprachen integriert werden.

Die Auswahl eines geeigneten Lösers hängt von den Systemanforderungen und der Systemeinrichtung ab. Die Löser mit fester Schrittweite sind sehr einfach zu implementieren und nicht an eine spezifische Sprache oder ein bestimmtes Tool gebunden, wohingegen die Löser mit variabler Schrittweite an eine Sprache bzw. ein Tool gebunden sind oder als Bibliothek integriert werden müssen – dafür sind sie aber stabiler und schneller.

Diese Methoden ermöglichen IMT, die Umsetzbarkeit komplexer thermodynamischer Systeme zu analysieren und mögliche Konstruktionsprobleme zu identifizieren. Darüber hinaus ermöglichen sie es IMT, eine Reglerkonstruktion effizient zu erarbeiten, selbst wenn ein Prototyp noch nicht zur Verfügung steht.