Matroids Matheplanet Forum Index
Moderiert von Wally haerter
Differentialgleichungen » Partielle DGL » Druckgradienten in einem Gasnetz
Autor
Beruf J Druckgradienten in einem Gasnetz
Jqb8314
Junior Letzter Besuch: im letzten Quartal
Dabei seit: 29.03.2018
Mitteilungen: 9
  Themenstart: 2022-06-06

Liebe Bewohner des Matheplaneten, ich möchte berechnen, wie sich der Druck in einem Gasrohr ändert, durch das Gas turbulent bei geringen Drücken fließt. Anders als beim klassischen Fourier/Fick'schen Gesetz ist der Druckverlust proportional zum Quadrat des Flusses. Ich komme auf folgende PDE (pdiff(u,t))^2 =-a*u(t,x)*(\partial^3 u)/(\partial x^3) Die Variable u ist die Gas Konzentration in mol/m3, also ein Maß für die Dichte oder den Druck im Gas. a ist eine positive Konstante. Die dritte Ableitung von u nach x sollte negativ sein. Das Ganze muss dann als eine Randwert-Aufgabe gelöst werden (Fluss am Ein- und Ausgang des Rohres ist bekannt, Druck-profil bei t = 0 ist bekannt). Sieht erstmal nicht schwierig aus. Wie löst man eine solche PDE? Ich habe im Studium leider nie PDEs behandelt. Sie sieht nicht aus wie eine Standard Elliptische, Parabolische oder Hyperbolische PDE, für die man Lösungen findet. Gibt es dafür überhaupt eine analytische Lösung? Viele Grüße JQB


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.1, eingetragen 2022-06-06

Hallo Jqb8314, \quoteon(2022-06-06 12:19 - Jqb8314 im Themenstart) Sieht erstmal nicht schwierig aus. Wie löst man eine solche PDE? Ich habe im Studium leider nie PDEs behandelt. \quoteoff Der Satz des Jahres. 😂 Also EINE Lösung zu finden ist hier nicht so schwer. ALLE Lösungen zu finden ist erheblich schwieriger. Diese DGL ist nichtlinear, wenn Du sie denn richtig hergeleitet hast. Das bedeutet, Superposition klappt nicht. Man könnte hier einen Produktansatz versuchen, und zwar $$u(t,x)=f(x)g(t)$$Dann folgt aus der PDE: $$f(x)\ddot g(t)=-af(x)g(t)\;f'''(x)g(t)$$$$\frac{\ddot g(t)}{g(t)^2}=-af'''(x)=c_1$$Demnach ist $$f(x)=-\frac{c_1}{6a}x^3+c_2x^2+c_3x+c_4$$Für $g(t)$ ergibt sich die DGL: $$\ddot g(t)=c_1g(t)^2$$Mit $\dot g(t)$ multiplizieren und integrieren: $$\dot g(t)\ddot g(t)=c_1\dot g(t)g(t)^2$$$$\tfrac12\dot g(t)^2=\tfrac13c_1g(t)^3+c_5$$$$\dot g(t)^2=\tfrac23c_1g(t)^3+2c_5$$Die Lösung dieser DGL sind die Weierstrass-P-Funktionen: $$g(t)=\wp\left(\sqrt{\frac{c_1}6}\;t;0,-12\frac{c_5}{c_1}\right)$$Dabei sind die Werte $g_2=0$ und $g_3=-12\frac{c_5}{c_1}$ Invarianten, die sich aus den Gitterkonstanten der doppeltperiodischen Weierstrass-P-Funktion ergeben. Die Umrechnung ist im Allgemeinen nicht ganz einfach, man muss zunächst die Nullstellen des Polynoms $4z^3-g_2z-g_3=0$ bestimmen, und aus den drei Nullstellen mittels arithmetisch-geometrischen Mittels die Gitterkonstanten berechnen. Die meisten CAS können das. Wenn Du nicht über die Weierstrass-P-Funktionen Bescheid weißt, wird das hier echt anstrengend werden... Aber diese Lösung ist nur eine spezielle Lösung, nicht die allgemeine. Wenn obiges $f(x)$ nicht genügt, um die Anfangsbedingung zu erfüllen, sehe ich wenig Chancen, eine allgemeine Lösung zu finden, wenn diese spezielle schon so kompliziert ist. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.2, eingetragen 2022-06-06

... wobei ich beim nochmal nachlesen mich gerade frage, ob es auf der linken Seite $$\left(\frac{\partial u}{\partial t}\right)^2$$(wie es da steht), oder $$\left(\frac{\partial}{\partial t}\right)^2 u$$heißen soll (wie ich gerechnet habe). Sollte die PDE stimmen, wie sie oben steht, was ich jetzt mal unterstelle, dann wäre meine Lösung falsch. Ciao, Thomas


   Profil
Jqb8314
Junior Letzter Besuch: im letzten Quartal
Dabei seit: 29.03.2018
Mitteilungen: 9
  Beitrag No.3, vom Themenstarter, eingetragen 2022-06-07

Vielen Dank für die ausführliche Antwort. Normalerweise gibt es auf dem Matheplaneten immer nur kryptische Hinweise und der Rest ist dann trivial (für den Antwortenden). Wenn man immer nur Angst vor komplizierten Formeln hat, kommt man nicht weit. 😉 Es ist nun mal so, dass Probleme im echten Leben meist keine einfache Lösung haben. Die Herleitung der Formel ist wie folgt: der Druckverlust einer stationären Strömung ist dp/dx =-(\rho v^2)/2 \lambda/D aus dem Idealen Gasgesetz und Bilanz um ein Volumenelement folgt pV = nRT, M = m/n, \rho = m/V, V = \Delta x A \rho = (nM)/(\Delta x A) v = n^*(\Delta x)/n dp = dn RT/(A\Delta x) Alles zusammen dn/dx =- \lambda/(2D) (M\Delta x^2)/(nRT) n^*^2 Ähnlich wie beim Fick'schen Gesetz definieren wir Flux und Konzentration u = n/(A \Delta X) J = n^*/A ODE mit Flux und Konzentration du/dx =-\lambda/(2D) M/uRT J^2 J = sqrt(-(2DuRT)/(\lambda M) du/dx) Massenerhaltung (\partial u)/(\partial t) +\partial/(\partial x) J = 0 Nun kommt das worauf Du hinweist (\partial u)/(\partial t) + \partial/(\partial x) sqrt(-(2DuRT)/(\lambda M) (\partial u)/(\partial x)) = 0 -(\partial u)/(\partial t) = \partial/(\partial x) sqrt(-(2DuRT)/(\lambda M) (\partial u)/(\partial x)) ((\partial u)/(\partial t))^2 = -au (\partial^3 u)/(\partial x^3) Stimmt die Herleitung? mit p [Pa] Druck x [m] Wegkoordinate \rho [ kg/m^3 ] Dichte v [m/s] Geschwindigkeit \lambda [-] Rohrreibungszahl (Darcy) V [m^3] Volumen n [mol] Stoffmenge R [ J/(mol K) ] allgemeine Gaskonstante T [K] Temperatur M [ kg/mol ] Molmasse m [kg] Masse \Delta x [m] Rohrabschnitt A [m^2] Rohrquerschnitt n^* [ mol/s ] Stoffmengenstrom u [ mol/m^3 ] Stoffmengendichte J [ mol/m ] Stoffmengenstromdichte a [ m^3/s^2 ?] Konstante (positiv) Die Lösung sollte real sein und monoton fallend in x. Das hat mir ein bisschen Hoffnung gemacht, dass die PDE loesbar sein koennte.


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.4, eingetragen 2022-06-07

Hallo Jqb8314, nein, die Herleitung ist wohl auf mehreren Ebenen falsch. Zunächst mathematisch an dieser Stelle: -(\partial u)/(\partial t) = \partial/(\partial x) sqrt(-(2DuRT)/(\lambda M) (\partial u)/(\partial x)) Hier kannst Du nicht einfach quadrieren, indem Du den $\frac{\partial}{\partial x}$-Operator für sich genommen quadrierst. Wenn der Strich für Ableitung nach $x$ stehe, dann ist $$\left(\frac{\partial}{\partial x}f\right)^2=f'^2\\ \left(\frac{\partial}{\partial x}\right)^2(f^2)=(f^2)''=(2ff')'=2f'^2+2ff''$$offensichtlich nicht das gleiche. Richtig wäre in der DGL oben, erst die Wurzel abzuleiten, wodurch selbige in den Nenner rutscht, und darüber hinaus müsstest Du das Innere der Wurzel nach der Produktregel ableiten (denn darin steht ja eigentlich $u\cdot u'$, warum solltest Du nur $u'$ ableiten?), und DANN könntest Du quadrieren. Dadurch würde die ganze DGL erheblich komplizierter. Aber auch an der hier zitierten DGL habe ich so meine Zweifel. Du möchtest eine Rohrströmung berechnen, aber was für eine Rohrströmung soll dass sein? Gut, turbulent ist klar, aber soll das ganze stationär ablaufen, d.h. Druck am Eintritt und am Austritt sind bekannt, und dann stellt sich ein gewisser stationärer Geschwindigkeitsverlauf ein? So könnte man ja vielleicht rechnen, aber dann hätte man gar keine Zeitabhängigkeit, der Verlauf des Drucks und der Geschwindigkeit entlang des Rohres wäre zeitlich konstant. Außerdem hat Gas nun einmal die Eigenschaft, dass seine Temperatur auch vom Druck abhängt, wenn eine recht plötzliche Druckänderung eintritt. Wie soll damit umgegangen werden? Strömt das Gas adiabat, also ohne Austausch (Abgabe oder Aufnahme) von Energie mit der Rohrwand? Oder ist die Strömung isotherm? Davon abhängig müsste man die Beziehung $pV^n=\text{konstant}$ mit ins Spiel bringen. Isothermische Strömung halte ich aber für unwahrscheinlich, da es eine vollständige Temperaturangleichung an die Rohrwand voraussetzt, was bei Gas eigentlich nur bei sehr langsamer Strömung der Fall wäre. Last but not least ist außerdem die Rohrreibungszahl $\lambda$ keine Konstante, sondern abhängig von der Reynoldszahl, die selbst geschwindigkeits-, druck- und viskositätsabhängig ist, und die Viskosität ist temperaturabhängig. $\lambda$ kannst Du also gar nicht als Konstante betrachten, und $T$ sehr wahrscheinlich auch nicht. Erzähl doch vielleicht ein bisschen über den Hintergrund der Berechnung, denn ich fürchte, wir müssen hier ganz vorne anfangen. Ciao, Thomas


   Profil
Jqb8314
Junior Letzter Besuch: im letzten Quartal
Dabei seit: 29.03.2018
Mitteilungen: 9
  Beitrag No.5, vom Themenstarter, eingetragen 2022-06-07

Hallo Thomas, Danke für die Korrektur. Ich werde die Produktregel das naechste Mal nicht vergessen. Zur eigentlichen Aufgabenstellung: Ich möchte den instationären Druckverlauf in einem langen Rohr vorhersagen, das an einem Ende mit Gas gefüllt wird und an dessen anderem Ende Gas entnommen wird. Gasbefüllung und Gasentnahme müssen nicht gleich sein. Es geht darum, einen Regelalgorithmus zu entwickeln, um die Gasentnahme nach dem Druck am Ausgang der Rohrleitung zu regeln. Dazu muss die Dynamik des Systems untersucht werden. Insbesondere der Anfahrvorgang (Rohr hat überall den gleichen Druck bei t = 0) ist von Bedeutung: Wie schnell steigt der Druck am Ende des Rohres, wenn dort keine Entnahme stattfindet aber am anderen Ende des Rohres Gas zugeführt wird? Ich mache ein paar Annahmen. 1. Die Strömung ist isotherm. Es handelt sich um ein langes Niederdruck-Gassnetz. Die Leitung ist in der Erde und nicht isoliert. 2. Die Rohrreibungszahl ist zwar abhängig von der Reynolds-Zahl, die Änderung ist aber klein. Man kann eine mittlere Rohrreibungszahl berechnen. 3. Die Druckschwankungen sind klein, der Druck wird zwischen 2 und 1.5 bar schwanken. Die Strömung kann als inkompressibel angenommen werden. Annahmen 1 und 2 sind meines Erachtens ok, bei 3 bin ich mir nicht so sicher. Ich werde es mal mit einer inkompressiblen Strömung versuchen. Vielleicht kürzen sich dann ein paar Terme weg und die Gleichung wird einfacher (Die Hoffnung stirbt zuletzt). Wenn das nicht funktioniert, kann ich das Ganze diskretisieren und die Leitung in 100 kleine Rührkessel aufteilen. Dann macht man aus der PDE 100 ODE – obwohl mich das Lösen einer PDE schon gereizt hätte. Viele Grüße Jürgen


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.6, eingetragen 2022-06-07

Hallo Jqb8314, EDIT: mit Deinem letzten Beitrag ist die nachfolgende Rechnung hinfällig. ich habe das mal durchgerechnet, unter der nassforschen Annahme, dass $\lambda$ als konstant angenommen werden kann, weil der Druckverlust hoffentlich nicht all zu hoch ist und die Änderungen am Zustand des Gases daher nicht gravierend sind. Normalerweise rechnen wir keine kompletten Lösungen vor. Das hier ist eine Ausnahme, weil es sich um eine Frage aus dem beruflichen Umfeld handelt und keine Übungsaufgabe für die Uni ist. Wir nehmen mal an, der Rohrquerschnitt sei konstant $A$, der Massenstrom sei konstant und zwar $\dot m$. Dann gilt aufgrund der Kontinuitätsgleichung: $$\varrho A v=\dot m$$$$\varrho v=\frac{\dot m}{A}\tag1$$Außerdem gilt $$p'=-\frac{\lambda}{2D}\varrho v^2\tag2$$Wir nehmen an, dass sich das Gas adiabat durch das Rohr bewegt, so dass also zu jeder Zeit $$\frac{p_0}{\varrho_0^\kappa}=\frac{p}{\varrho^\kappa}\tag3$$$\kappa$ ist der konstante, aber gasabhängige Adiabatenexponent (für zweiatomige Gase gilt $\kappa=1\mathord,4$), die Größen $p,\varrho$ und $v$ seien jeweils abhängig von der Ortskoordinate $x$, aber nicht abhängig von $t$ (stationäre Strömung), der Strich bedeute Ableitung nach $x$, der Index 0 beziehe sich auf den Zustand am Eintritt in das Rohr. Mehr als das hier brauchst Du nicht. Du kannst es nun selbst versuchen, auszurechnen, oder Du kannst nachfolgend "spicken". \showon Nun muss man nur ein bisschen ineinander einsetzen. Zunächst multiplizieren wir (2) mit $\varrho$: $$p'\varrho=-\frac{\lambda}{2D}\varrho^2 v^2\tag4$$und setzen hier rechts die Gleichung (1) ein: $$p'\varrho=-\frac{\lambda\dot m^2}{2DA^2}\tag5$$(3) lösen wir nach $\varrho$ auf: $$\varrho=\varrho_0p_0^{-\frac1\kappa}\;p^{\frac1\kappa}\tag6$$und setzen in (5) ein: $$\varrho_0p_0^{-\frac1\kappa}\;p'p^{\frac1\kappa}=-\frac{\lambda\dot m^2}{2DA^2}\tag7$$Bissi aufräumen: $$p'p^{\frac1\kappa}=-\frac{\lambda\dot m^2 p_0^{\frac1\kappa}}{2DA^2\varrho_0}\tag8$$Rechts stehen nur Konstanten (zumindest tun wir mal so), daher können wir nun integrieren: $$\frac{p^{\frac1\kappa+1}}{\frac1\kappa+1}=c_1-\frac{\lambda\dot m^2 p_0^{\frac1\kappa}}{2DA^2\varrho_0}x\tag9$$Es muss gelten $p(x=0)=p_0$, so dass wir für die Integrationskonstante direkt einsetzen können: $$\frac{p^{\frac1\kappa+1}}{\frac1\kappa+1}=\frac{p_0^{\frac1\kappa+1}}{\frac1\kappa+1}-\frac{\lambda\dot m^2 p_0^{\frac1\kappa}}{2DA^2\varrho_0}x\tag{10}$$$$p^{\frac1\kappa+1}=p_0^{\frac1\kappa+1}-\frac{\lambda\dot m^2 p_0^{\frac1\kappa}(\frac1\kappa+1)}{2DA^2\varrho_0}x\tag{11}$$Und damit: $$p(x)=\left(p_0^{\frac{\kappa+1}\kappa}-\frac{\lambda\dot m^2 p_0^{\frac1\kappa}(\kappa+1)}{2\kappa DA^2\varrho_0}x\right)^{\frac{\kappa}{\kappa+1}}\tag{12}$$Noch ein bisschen schöner: $$p(x)=p_0\left(1-\frac{(\kappa+1)\lambda\dot m^2}{2\kappa DA^2\varrho_0p_0}x\right)^{\frac{\kappa}{\kappa+1}}\tag{13}$$ Versuche, das mal nachzuvollziehen. Solltest Du Fragen zu den einzelnen Gleichungen haben, kannst Du gerne nachfragen. \showoff Ciao, Thomas [Die Antwort wurde nach Beitrag No.4 begonnen.]


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.7, eingetragen 2022-06-07

Hallo Jürgen, \quoteon(2022-06-07 13:37 - Jqb8314 in Beitrag No. 5) (...)Insbesondere der Anfahrvorgang (Rohr hat überall den gleichen Druck bei t = 0) ist von Bedeutung: Wie schnell steigt der Druck am Ende des Rohres, wenn dort keine Entnahme stattfindet aber am anderen Ende des Rohres Gas zugeführt wird? (...) Die Strömung kann als inkompressibel angenommen werden. \quoteoff Das passt nicht zueinander. Wenn das Medium inkompressibel ist, dann ist die Druckerhöhung auch "sofort" da., und dann kannst Du nicht an einem Ende einfüllen, während Du am anderen Ende nichts entnimmst. Das geht nur mit kompressiblem Fluid. Wenn Du in einer Gasleitung an einem Ende schlagartig den Druck erhöhst, indem Du z.B. ein Ventil öffnest, dann breitet sich die Druckerhöhung mit Schallgeschwindigkeit aus, und es wird sicher zu Wellenphänomenen kommen. Man kann zwar die Wellengleichung für den Druck in dem Rohr berechnen, aber dabei wird von einem stillstehenden Medium ausgegangen, nicht von einem strömenden. Das, was Du berechnen willst, ist daher schon tricky. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.8, eingetragen 2022-06-08

Hallo Jürgen, wenn wir das Ganze instationär betrachten wollen, müssen wir meine Gleichungen oben etwas erweitern: Nach wie vor gilt adiabatische Umwandlung. Das heißt, wir nehmen an, dass das Gas in einem bestimmten Zustand in die Leitung gepumpt wird und diesen Zustand energetisch beibehält, also keine Zeit hat, Energie mit der Rohrwand auszutauschen. Einen Tod kann man nur sterben. Wenn man $pV^n=\text{const}$ ansetzt, dann ist $n=\kappa$ bei adiabatischen Zustandsänderungen, und $n=1$ bei isothermen Zustandsänderungen. Eventuell kann man einen Polytropenexponenten von $n=1\mathord,2$ annehmen, um einen limitierten Energieaustausch zu improvisieren. Aber wie gesagt, wir bleiben erst einmal adiabat. Das heißt, es gilt nach wie vor: $$\frac{p_0}{\varrho_0^\kappa}=\frac{p}{\varrho^\kappa}\tag1$$Darüber hinaus gilt wegen der Quellenfreiheit: $$\frac{\partial \varrho}{\partial t}+\frac{\partial}{\partial x}(\varrho v)=0\tag2$$(Wenn man hier im stationären Fall die Zeitableitung null setzte, würde daraus die Kontinuitätsgleichung folgen). Außerdem benötigen wir noch den Energieverlust aufgrund der Rohrreibung. Hier hast Du allerdings die Gleichung für den stationären Fall verwendet, was bedeutet, dass Du unterschlagen hast, dass eine Druckdifferenz, die größer ist als der Reibungsverlust, zu einer Beschleunigung des Gases führt. Das bedeutet: $$\frac{\partial p}{\partial x}+\frac{\lambda}{2D}\varrho v^2+\varrho\frac{\partial v}{\partial t}=0\tag3$$Mit Gleichung (1) kann man $p$ noch verhältnismäßig leicht aus Gleichung (3) eliminieren: $$\frac{p_0}{\varrho_0^\kappa}\frac{\partial}{\partial x}(\varrho^\kappa)+\frac{\lambda}{2D}\varrho v^2+\varrho\frac{\partial v}{\partial t}=0\tag4$$Die beiden partiellen Differentialgleichungen (2) und (4) beschreiben nun weitestgehend den instationären Fall. Ich sehe es allerdings im Moment als aussichtslos an, nachdem ich gestern Abend schon ein Weilchen darauf herumgekaut habe, diese Gleichungen irgendwie sinnvoll ineinander einzusetzen und eine lösbare partielle Differentialgleichung entweder für $\varrho$ oder für $v$ zu erhalten. Vielleicht haben Mathematica oder Maple dazu eine Idee. Ciao, Thomas


   Profil
Jqb8314
Junior Letzter Besuch: im letzten Quartal
Dabei seit: 29.03.2018
Mitteilungen: 9
  Beitrag No.9, vom Themenstarter, eingetragen 2022-06-08

Hallo Thomas, zuerst einmal vielen Dank für Deine Zeit. Ich weiß das sehr zu schätzen. Die Sache mit der Kompressibilität stimmt leider. Dein Ansatz unten mit der erweiterten Darcy Gleichung für eine beschleunigte Strömung (Gleichung 3) ist sehr hilfreich und war bisher immer eine falsche Annahme von mir. Ich habe Deine Gleichung 4 für den isothermen Fall vereinfacht und versucht zu lösen, aber man kann leider nicht eine der beiden Variablen eliminieren, sondern muss Fluss und Konzentration (oder Geschwindigkeit und Dichte) simultan lösen. Ich werde es mal numerisch versuchen und Matlabs PDE-solver kaufen. Ich habe schon mal mit Randwert-solvern in MATLAB gearbeitet und traue mir diesen Lösungsweg eher zu. Mein Mathelehrer hat immer gesagt: Mathematik ist durch Denken Rechnen zu vermeiden, aber manchmal kommt man mit einem numerischen Ansatz schneller weiter. Am Ende zählt das Ergebnis. Vielen Dank für Deine Hilfe Juergen


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.10, eingetragen 2022-06-08

Hallo Jürgen, ich könnte mir sogar vorstellen, dass man im ersten Anlauf die Rohrreibung vernachlässigen kann. Die wird am Ende zwar dafür sorgen, dass sich überhaupt ein stationärer Zustand einstellt, aber für die Frage, wann die Druckerhöhung am Leitungsende ankommt, kann man sie vielleicht vernachlässigen. Wenn man außerdem isotherm rechnet, was natürlich auch eine verwegene Vereinfachung ist, dann käme man auf folgendes, etwas einfacheres System von Differentialgleichungen: $$\frac{\partial \varrho}{\partial t}+\frac{\partial}{\partial x}(\varrho v)=0\tag2$$und $$\frac{p_0}{\varrho_0}\cdot\frac{\partial \varrho}{\partial x}+\varrho\frac{\partial v}{\partial t}=0\tag{4b}$$Definieren wir nun eine Größe $f$, für die ich nicht einmal einen sinnvollen physikalischen Namen habe, mit $$\frac{\partial f}{\partial x}=\varrho$$und setzen das in (2) ein, dann folgt: $$\frac{\partial^2 f}{\partial t\partial x}+\frac{\partial}{\partial x}(\varrho v)=0\tag5$$Diese Gleichung können wir nach $x$ integrieren: $$\frac{\partial f}{\partial t}+\varrho v=c(t)\tag6$$Dabei ist $c(t)$ eine Integrationskonstante (in $x$), die durchaus von $t$ abhängen kann, aber die Ableitung nach $x$ wäre halt null, damit (5) erfüllt ist. Diese Funktion müsste man später durch Anfangs- oder Randbedingungen bestimmen. Leiten wir diese Gleichung nach $t$ ab: $$\frac{\partial^2 f}{\partial t^2}+v\frac{\partial \varrho}{\partial t}+\varrho\frac{\partial v}{\partial t}=\dot c(t)\tag7$$Davon ziehen wir (4b) ab: $$\frac{\partial^2 f}{\partial t^2}-\frac{p_0}{\varrho_0}\cdot\frac{\partial \varrho}{\partial x}+v\frac{\partial \varrho}{\partial t}=\dot c(t)\tag8$$Nun multiplizieren wir diese Gleichung mit $\varrho$, Gleichung (6) mit $\frac{\partial \varrho}{\partial t}$, und ziehen das voneinander ab, um $v$ zu eliminieren: $$\frac{\partial^2 f}{\partial t^2}\varrho-\frac{p_0}{\varrho_0}\cdot\frac{\partial \varrho}{\partial x}\varrho-\frac{\partial f}{\partial t}\frac{\partial \varrho}{\partial t}=\dot c(t)\varrho-c(t)\frac{\partial \varrho}{\partial t}\tag9$$Wir ersetzen nun überall noch $\varrho$, und erhalten: $$\frac{\partial^2 f}{\partial t^2}\cdot\frac{\partial f}{\partial x}-\frac{p_0}{\varrho_0}\cdot\frac{\partial^2 f}{\partial x^2}\cdot\frac{\partial f}{\partial x}-\frac{\partial f}{\partial t}\frac{\partial^2 f}{\partial t\partial x}=\dot c(t)\frac{\partial f}{\partial x}-c(t)\frac{\partial^2 f}{\partial t\partial x}\tag{10}$$Das ist ein kleines bisschen unübersichtlich. Verwenden wir den Punkt für die Ableitung nach $t$ und den Strich für die Ableitung nach $x$, dann haben wir: $$\left(\ddot f-\frac{p_0}{\varrho_0}f''\right)f'-\dot f\dot f'=\dot c(t)f'-c(t)\dot f'\tag{11}$$Macht's auch nicht besser. Wir haben nun eine einzige partielle Differentialgleichung für $f$, woraus wir mittels Ableitung nach $x$ anschließend die Dichteverteilung $\varrho$ bestimmen könnten, müssten aber auch noch irgendwie die Integrationsfunktion $c(t)$ bestimmen, die zeit- aber nicht ortsabhängig ist. Und durch die Tatsache, dass $f$ links immer in zweiter Potenz auftaucht, ist die DGL auch noch nichtlinear. Man erkennt allerdings in der großen Klammer einen Term, der typisch ist für eine Wellengleichung, womit schon fest steht, dass das Füllen der Rohrleitung mit Gas nicht schön gleichmäßig und monoton (also Druck und Dichte an jedem Punkt streng monoton steigend mit der Zeit), sondern wellenförmig passiert, also "schwappend". War ja irgendwie zu erwarten. Eine analytische Lösung ist wohl wirklich nicht drin, selbst für diesen vereinfachten Ansatz nicht. Ich denke, an einer numerischen Lösung führt kein Weg vorbei. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.11, eingetragen 2022-06-11

Hallo Jürgen, ich habe mittels der Wolfram Cloud mal eine numerische Lösung berechnet. Ergebnisse siehe unten. Zunächst habe ich die Gleichungen in eine dimensionslose Form überführt. Mit $c^2=\kappa\frac{p_0}{\varrho_0}$ Schallgeschwindigkeit, $l,d$ Länge und Durchmesser des Rohres, $\lambda$ Rohrreibungszahl $k=\frac{\lambda l}{2d}$ dimensionslose Widerstandszahl, $\xi=\frac xl$ dimensionslose Rohrkoordinate, $T=\frac lc$ Laufzeit des Schalls durch das Rohr, $\tau=\frac tT$ dimensionslose Zeit, $r=\frac{\varrho}{\varrho_0}$ relative Dichte ($r=1$ entspricht der Dichte im Rohr bei $t=0$), $m=\frac vc$ relative Geschwindigkeit (Mach-Zahl) erhält man folgende DGL: $$\frac{\partial r}{\partial \tau}+\frac{\partial}{\partial \xi}(r\,m)=0$$und $$r^{\kappa-2}\frac{\partial r}{\partial \xi}+k\,m|m|+\frac{\partial m}{\partial \tau}=0$$Bei der letzten PDGL muss man beim Strömungswiderstand statt $v^2$ mit $v|v|$ rechnen (hier $m|m|$), um der Richtungsumkehr der Reibung mit der Geschwindigkeit Rechnung zu tragen. Wie man sieht, hängt die Lösung der PDGL nun nur von dem Adiabatenexponenten $\kappa$ und der dimensionslosen Widerstandszahl $k$ ab. Bei einer Leitungslänge von 100m, Durchmesser 1cm und einer Reibungszahl von $\lambda=0\mathord,02$ gilt zum Beispiel $k=100$, was die zweite DGL ziemlich dominiert. Nachfolgend die Lösungen für $\kappa=1\mathord,4$ (zweiatomiges Gas), $k=100$, Dichteerhöhung um den Faktor $1\mathord,2$ (z.B. von 2,5 auf 3 bar bei konstanter Temperatur). Dichteverlauf: Im Vordergrund ist die Zeitachse, im Hintergrund die Raumachse. Man erkennt schön folgende vorhersehbare Effekte: a) Die Dichteerhöhung (und damit auch die Druckerhöhung) kommt erst nach $\tau=1$ am Rohrende an, die Diagonale von "links" nach "rechts" ist die Wellenfront, die sich durch das Rohr arbeitet. b) Kurz hinter dem Eintritt in das Rohr ("hinten" bei $\xi\approx 0$) steigt die Dichte in einer rasch abklingenden Schwingung auf den höheren Wert, der bei $\xi=0$ als Randbedingung vorgegeben wird. Durch die starke Dämpfung gibt es keine sehr starken Welleneffekte. Geschwindigkeitsverlauf: Hier erkennt man noch deutlich die recht stark schwingende Zuflussgeschwindigkeit am Rohreintritt. Am fernen (und verschlossenen) Rohrende im Vordergrund rechts muss die Geschwindigkeit null sein (Randbedingung). Wenn die Wellenfront am Rohrende ankommt, liegt über die komplette Rohrlänge eine recht gleichmäßige Geschwindigkeitsverteilung vor (linker Vordergrund). Betrachten wir nun einen längeren Zeitraum, nämlich bis $\tau = 12$ (was bei Luft in einem 100m langen Rohr knappe 4 Sekunden bedeutet). Zunächst wieder der Dichteverlauf: Nachdem der recht heftige Druckausgleich zu Beginn (rechter Rand) und die Ausbreitung mit Schallgeschwindigkeit abgeschlossen sind, stellt sich eine relativ langsame Schwingbewegung ein. Im ganzen Rohr schwingt die Dichte mit einer Periode von etwa $4T$ (rund 1,3 Sekunden bei eben genannten Beispielrohr). Geschwindigkeitsverlauf: Durch die langsame Schwingung fällt die Geschwindigkeit gelegentlich knapp unter null, so dass also Gas nicht nur zuströmt, sondern in geringen Mengen immer wieder auch zurückströmt, aber nach einer Dauer von rund $6T$ oder in diesem Beispiel mit 100m Rohrlänge nach rund zwei Sekunden hat sich die Gasströmung weitestgehend beruhigt. Ciao, Thomas


   Profil
MontyPythagoras
Senior Letzter Besuch: in der letzten Woche
Dabei seit: 13.05.2014
Mitteilungen: 2954
Wohnort: Werne
  Beitrag No.12, eingetragen 2022-06-12

Hallo nochmal, hier noch eine Animation der Dichte und der Geschwindigkeit bis $\tau=6$, um das Ganze noch ein wenig anschaulicher zu machen: Gasdichte: https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_GasDichte.gif Man erkennt schön, wie das Gas in das Rohr hineinschwappt, es gibt auch eine rücklaufende Welle vom verschlossenen Rohrende. Gasgeschwindigkeit: https://matheplanet.de/matheplanet/nuke/html/uploads/b/39826_GasGeschw.gif Zu Beginn findet der recht heftige Druckausgleich statt, was zu einer kurzfristigen Gasschwingung führt. Interessant finde ich, dass vor der Druckwelle ein Wellental herläuft, wo die Geschwindigkeit knapp negativ ist - wie bei einem Tsunami, wo das Wasser am Strand erst abfließt, bevor die Welle kommt. Ich bin nicht sicher, ob das real ist, oder ein Artefakt der numerischen Berechnung, weil die Randbedingungen recht krass wechseln und es zumindest in den Ableitungen der Randbedingungen Unstetigkeiten gibt. Ciao, Thomas


   Profil
Jqb8314
Junior Letzter Besuch: im letzten Quartal
Dabei seit: 29.03.2018
Mitteilungen: 9
  Beitrag No.13, vom Themenstarter, eingetragen 2022-06-12

Hallo Thomas, ich sehe, Dich lässt das Problem auch nicht los. Ich habe mal ein bisschen Literatur gewälzt. Es gibt mehrere Lösungswege. Bei einer isothermen Strömung gibt es eine analytische Lösung. (Gugat und Wintergerst) Dabei wird als Bezugssystem eine konstante Geschwindigkeit vorgegeben (normalerweise Schallgeschwindigkeit, die ist hier aber nicht konstant, da das Realgasverhalten berücksichtigt wird.) Die Lösung benötigt die Lambert W-Funktion auf der negativen Achse. Ist also nur semi-analytisch. Viele andere Literaturstellen lösen das nicht-isotherme Problem mit Bilanzen für Masse, Impuls und Energie, z.B. Sanaye und Mahmoudimehr, Ferreira et al, Sund und noch einige andere. Die zeitlichen und örtlichen Ableitungen werden als finite Differenzen berechnet, entweder mit Backward Time, Centered Space Method oder einer Mischung aus BTCS und Crank–Nicolson. Ich habe diesen Ansatz in Matlab umgesetzt und dabei alle Ableitungen händisch programmiert. Man muss eine 3N,3N-Matrix invertieren, die aus 9 tridiagonalen Untermatrizen besteht, (N = Anzahl der diskreten Rohrleitungssegmente), das Ganze ist sehr viel Code und irgendwo ist noch ein Fehler. Dann habe ich Deine Lösung mit der Wolfram-Cloud gesehen und dachte, das muss auch einfacher gehen. Man kann die Gleichungen umstellen, so dass man drei Gleichungen hat, bei denen jeweils eine zeitliche Ableitung auf der linken Seite steht und die örtlichen Ableitungen und der Rest auf der rechten Seite. Diese Form kann man in Matlabs PDE-solver übergeben, aber leider braucht der Solver spezielle Form der Randwerte. Das Problem hat zwei Randwerte im Fluss (Ein- und Aufgang) und einen in der Temperatur und das ist nicht kompatibel mit der Struktur des Matlab solvers. Im Paper von Chaczykowski werden die drei Gleichungen so umgestellt, dass sich drei Formeln für drei zeitliche Ableitungen ergeben. Es werden nur die örtlichen Ableitungen mit finiten Differenzen berechnet und die zeitlichen Ableitungen in ODE überführt. Das Ganze kann dann mit Solvern für steife ODE gelöst werden (so der Autor). Ich werde das heute Abend mal versuchen. Die Animationen sind sehr hilfreich, um den Regler auszulegen. Man muss sich vermutlich Gedanken zur Eigenfrequenz der Leitung machen und ob es zu Resonanzen kommen kann wenn man die Entnahme am Ende nach dem Druck am Ende regelt. Viele Grüße Juergen


   Profil
Jqb8314 hat die Antworten auf ihre/seine Frage gesehen.
Jqb8314 hat selbst das Ok-Häkchen gesetzt.
Jqb8314 wird per Mail über neue Antworten informiert.

Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2022 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]