Numerik: besserer Ansatz für DGL-System

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Chris987

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 12.03.13

	Wohnort: ---

	Version: ---

Verfasst am: 23.04.2013, 14:59 Titel: Numerik: besserer Ansatz für DGL-System

Hallo!
In diesem Thema geht es um die Lösung eines (DGL)-Gleichungsystems in, wie ich finde, ungünstiger Form.
Ich möchte mich gerne Schritt für Schritt vorantasten bzw. optimieren. Deswegen möchte ich erstmal nur eine GL betrachten.

Die sieht wie folgt aus (aus irgendeinem grund kommt das nicht als math-Ausdruck):
$ \frac{dN}{dt} = c1 \cdot c2 \cdot c3 \cdot \sqrt{Druck} \cdot e^{-21100/Temperatur} \cdot c4} $

N bezeichnet eine Anzahl gebildeter Partikel. Die Gleichung beschreibt also gewissermaßen die zeitliche Bildungsrate dieser Partikel.
Für dieses Modell habe ich den Druck und Temperatur für einen bestimmten Zeitraum als diskrete Messwerte vorliegen. Bei den anderen Werten handelt es sich um Konstante Werte, die nicht weiter interessieren.

Für diesen Zeitraum, gewissermaßen die Integrationsgrenzen, möchte ich wissen, wie viele Partikel neu gebildet werden.

Ich habe aus der Gleichung mathematisch eine Differenzengleichung gemacht in dem Sinne:
$ \frac{dN}{dt}=..... \\ \frac{\Delta N}{\Delta t}=...\\ {\Delta N}=...\cdot {\Delta t} $

Meine erste Lösung des Problems war folgende:

Code:

P_neu= c1 * c2 * c3 * sqrt(Zylinderinnendruck{ink,1}) * exp(-21100/T_zyl{ink,1}) * c4 * sw_sekunden

Funktion ohne Link?

ink ist das Zeitinkrement, also der Wert des Drucks und der Temperatur zu einem bestimmten Zeitpunkt. Diese Werte werden in einem Zeitabstand "sw-sekunden" neu gemessen und berechnet. Das heißt für das Zeitintervall "sw_sekunden" nehm ich Druck und Temperatur als konstant und multipliziere das mit der Intervalldauer.
Ergebnis ist die Anzahl neu gebildeter Partikel in der Intervalldauer.

Nach entsprechender Kritik habe ich versucht, das ganze besser zu lösen.
Ich habe lange und ausgiebig nach Lösungsmöglichkeiten per dsolve und ode45 gesucht, aber nichts passendes für mein Problem gefunden.

Mein zweiter Ansatz sieht so aus:

Code:

P_neu= c1 * c2 * c3 .* sqrt(Zylinderinnendruck) .* exp(-21100/T_zyl) * c4

for w=1:(BD+1)

rp_neu_int(w,1) = trapz(x_zeit(w:w+1,1),rp_neu(w:w+1,1)); % pro zeitabschnitt gebildete partikel
end

Funktion ohne Link?

Diesmal keine Differenzengleichung mehr. Stattdessen benutze ich einen Stützwertvektor P_neu für die Gesamtzeitdauer BD mit den Stützstellen x_zeit (Schrittweite sw_sekunden) um die Pro Zeitabschnitt gebildeteten Partikel zu berechnen. Das Integral lasse ich mir mit trapz berechnen.

Nun meine Frage an euch: Wie kann ich das eleganter lösen?
Ich glaube, dass es mit den Bordmitteln von matlab gehen sollte. Habe aber leider nicht den Durchblick Sad

Falls ich etwas unverständlich formuliert habe sollte, bitte schreibt, damit die Problemstellung deutlicher formulieren kann.

Gruß
Christian

Bluesmaster

Forum-Century


	Beiträge: 203

	Anmeldedatum: 13.11.11

	Wohnort: Gera

	Version: 2012a

Verfasst am: 23.04.2013, 22:55 Titel:

Mein Vorschlag:

- Die Formel noch mal ordentlich (zur Not als Bild)

- Die Frage auf 50 Worte eindampfen

Ohne Experte zu sein glaube ich nicht dass es (einfache) DGLs
gibt die man mit ode nicht lösen könnte

Gruß

Blues

Matmarv

Forum-Fortgeschrittener


	Beiträge: 54

	Anmeldedatum: 12.03.13

	Wohnort: ---

	Version: Matlab R2013a

Verfasst am: 24.04.2013, 09:59 Titel:

Ich weiß nicht ob ich dich richtig verstanden habe.
Versucht du die DGL zu lösen?

Du könntest die DGL in Zustandraumdarstellung bringen. Dann kennst du die Änderung des Zustandes.

also

x= Zustandsvektor
A=Systemmatrix
B= Eingangsvektor
u=Eingangsgrößen
xpunkt=Ableitung des Zustandes
xalt = letzter Zustand k-1

xpunkt = A*x + B*u
x = xalt + xpunkt *deltat

Gruß
Marvin

Chris987

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 12.03.13

	Wohnort: ---

	Version: ---

Verfasst am: 25.04.2013, 11:58 Titel: Neuformulierung

Hallo,
danke für die Antworten. Ich werde die Gleichungen nochmal neu schreiben.

Partikelneubildung; Druck und Temperatur liegen als diskrete Werte über der Zeit vor (im Abstand von 0,00001s gibt es einen neuen Wert), der Rest ist konstant.
$\frac{dN}{dt}=c1\cdot c2\cdot \sqrt{Druck}\cdot e^{\frac{-21000}{Temperatur}} \cdot c4$

Partikeloxidation; Druck und Temperatur wie gehabt, Anteile von Oxidatoren(Konz) werden aus Messdaten berechnet und liegen ebenso als Vektor vor.
$\frac{dNox}{dt}=-k1 \cdot k2 \cdot Konz \sqrt{Temperatur}\cdot (\pi N)^{\frac{2}{3}} \cdot [\frac{6 Mpartikel}{k4}]*k5$
$\frac{dNo2}{dt}=-w1 \cdot e^{\frac{-19778}{Temperatur}} \sqrt{Temperatur} Konz2 \cdot (\pi N)^{\frac{2}{3}} \cdot [\frac{6 Mpartikel}{w4}]*w5$

1. Wie kann man die Bildung der Partikel geschickt(er) berechnen?
(bisherige Ansätze habe ich im Eröffnungspost beschrieben)

2. Partikel werden auch oxidiert. Dies ist Temperaturabhängig, konzentrationsabhängig und abhängig von der momentanen Partikelkonzentration. Die Partikelmasse wird aus der Anzahl berechnet (Mpartikel).

Bisher habe ich das System als Differenzengleichung behandelt.
Als Zeitintervall Delta-t habe ich die Zeitschrittweite von Temperatur/Druck genutzt. Für jeden Zeitschritt habe ich dann eine Bilanz gemacht, die Schleife dann so oft iteriert, wie der Prozess dauert:

Code:

R_partikel(1+ink,1) = rp_neu_int(ink+1,1) + N_oxO2_int(ink,1) + N_oxOH_int(ink,1)

Funktion ohne Link?

Wie kann man die Entwicklung der Partikelanzahl mit einem numerischen Werkzeug lösen?

Bluesmaster

Forum-Century


	Beiträge: 203

	Anmeldedatum: 13.11.11

	Wohnort: Gera

	Version: 2012a

Verfasst am: 25.04.2013, 15:31 Titel:

Ok soweit zum fachlichen ( bitte habe auch Verständniss, dass nicht jeder das System nachvollziehen kann/ will ).

Wir müssen weiter abstrahieren, um wirklich Matlab Hilfe leisten zu können.
Matlab solver können jede Art von Problem lösen d.h. wie Matmarv
gesagt hat: Aufschlüsseln abhängige, unabhängige Größen:

Ich würde sagen, wir brauchen mal einen der sich mit DGLs gut
auskennt um deine Größen in eine geeignete Form zu bringen.
Ich mache trotzdem mal den Anfang:

ODE-System der Kategorie "Randwertproblem" (?)

Y1' = f( p , T )
Y2' = f( k , T , Y1 )
Y2' = f( c , T , Y2 )

p = Druck
T = Temp
k = Konzentration Oxidatoren

...liegen als Zeitverlauf f(t) als Randwerte vor.

Aber das ist alles noch sehr unsicher. Was genau ist Nox?
Die Anzahl der oxidierten Partikel?
Was ist No2 die anzahl der Oxidatoren?

Und wo ist die Bilanz, dass ein oxidiertes Partikel von
der (Rein)Partikel anzahl abgezogen wird?

Hier ist noch viel zu tun um das in eine ordentliche Form zu bringen.
Wie gesagt wer Ahnung von DGLs hat bitte melden!!!

Gruß

Blues

Chris987

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 12.03.13

	Wohnort: ---

	Version: ---

Verfasst am: 25.04.2013, 18:40 Titel:

Mein Verdacht lag bisher auf einem Anfangswertproblem. Von DGL habe ich leider nicht übermäßig viel Ahnung.

Deine Schreibweise habe ich mal übernommen und abgeändert:

Y1' = f( p , T )
Y2' = f( Konz1 , T , Y1 )
Y3' = f( Konz2 , T , Y1 )

Y1' ist die positive Bildungsrate von Partikeln. Y2' und Y3' sind dann negative Bildungsraten (negatives Vorzeichen in den Formeln). Zur Unterscheidung heißen sie NoxO2 (Sauerstoff als Namensgeber) und NoxOH (Hydroxyl). Konz1 und 2 liegen genau wie p und T als Zeitverlauf f(t) vor.

Das Y1 kommt einmal als Anzahl N vor, und einmal als Masse (Mpartikel). Lässt sich aber leicht ineinander umwandeln durch Multiplikation mit einem konstanten Faktor.

Eine Bilanzgleichung ist in meiner Quelle gar nicht genannt.
Die hab ich mir dazu überlegen (müssen).

Die arbeitet als Schleife Abschnittsweise (keine exakte Lösung!!):

Zeitintervall1: Y1 = p_neu(1) - 0 - 0
Zeitintervall2: Y2 = Y1 + p_neu(2) -pox(Y1) - poh(Y1)
Zeitintervall3: Y3 = Y2 + p_neu(3) -pox(Y2) - poh(Y2)
Zeitintervall4: Y4 = ......

Danke schon einmal an alle, die sich hiermit beschäftigen und Gedanken machen!

Bluesmaster

Forum-Century


	Beiträge: 203

	Anmeldedatum: 13.11.11

	Wohnort: Gera

	Version: 2012a

Verfasst am: 25.04.2013, 21:03 Titel:

Hast recht, ist ein AWP,

mich haben die gegebenen Zeitverläufe irritiert, aber die Partikelzahl
selbst ist natürlich zu keinem Zeitpunkt bekannt.

Gute Nachricht:

Es ist gar kein DGL-System, sondern eine gewöhnliche DGL 1.Ordnung

Y' = f( Y , t ) der denkbar einfachste Fall.

Das liegt daran dass die Partikelsenken Nox und Noh nichts weiter sind
als negative Partikel also mathematisch die gleiche Variable.

Normalerweise mussen DLG-Solver ran, wenn:

- ich in meinem System nicht alles Messen kann was ich will
- Versuche zu aufwändig sind
- ich nicht alles alle Variablen so einstellen kann wie ich will

Hier fiele mir nur ein, dass du an den Konstanten rumfummeln kannst,
solange bis sie zur echten Partikelanzahl passen (die du doch garantiert auch messen kannst oder?)

Nun denn ich habe mal einen lauffähigen Rohling ausgebastelt,
läuft mit runge-kutta45. Versuche mal den nachzuvollziehen,
damit zu experiementieren und dich wieder zu melden.

Gruß

Blues

Code:

function numberOfParticles = particleSimulation( numberOfParticlesAtStartTime , endTime )
% calculates the creation of particles over time

% Testmode
if nargin < 1 , numberOfParticlesAtStartTime = 10e8 ; end % Falls keine Partikelanfangszahl gegeben
if nargin < 2 , endTime = 1 ; end % Falls keine Simulationszeit gegeben

nMeasurePoints = 10000000;

global measurementTimeStep measuredPressure measuredTemperature measuredconOh measuredconOx

measurementTimeStep = 0.00001; % [s]

% Prüfen ob die Endzeit im Messrahmen liegt
assert( endTime < measurementTimeStep * nMeasurePoints , 'given endTime exceeds the limits of measurements for pressure, temperature and concentration' );

% Messwerte vorübergehend auf konst = 1 setzen
measuredPressure = ones( nMeasurePoints , 1 );

measuredTemperature = ones( nMeasurePoints , 1 );

measuredconOh = ones( nMeasurePoints , 1 );

measuredconOx = ones( nMeasurePoints , 1 );

options = odeset( 'RelTol' , 1e-4 , 'AbsTol' , 1e-4 );
[ T , Y ] = ode45( @particleCreationRate ,[ 0 endTime ] , numberOfParticlesAtStartTime , options );

% [ T , Y ] = ode45( @particleCreationRate ,[ 0 endTime ] , numberOfParticlesAtStartTime )

numberOfParticles = Y(end);

% Testplot
plot(T, Y);

function dydx = particleCreationRate( t , y )

% Konstantendefinition: (übergangsweise alle = 1)
c1 = 1; c2 = 1; c3 = 1; c4 = 1; % c3 fehlt in Formel??
k1 = 1; k2 = 1; k3 = 1; k4 = 1; k5 = 1; % k3 fehlt in Formel??
w1 = 1; w2 = 1; w3 = 1; w4 = 1; w5 = 1; % w2/ w3 fehlen in Formel???
M = 10^-9; % [kg/Partikel]

dydx = + c1 * c2 * sqrt( pressure(t) ) * exp( -21000 / temperature(t) ) * c4 ... Neubildungsanteil
- k1 * k2 * conOx(t) * sqrt( temperature(t) ) * ( pi * y)^(2/3) * ( 6 * M * y / k4 ) * k5 ... Sauerstoffoxidationsanteil
- w1 * exp( -19778 / temperature(t) ) * sqrt( temperature(t) ) * conOh(t) * ( pi * y )^(2/3) * ( 6 * M * y / w4 ) * w5; % Hydroxyloxidationsanteil

function p = pressure( t )

global measurementTimeStep measuredPressure

% Interpolation möglich aber sinnlos bei so trägen Größen wie Temperatur, Druck und Konzentration
indexOfMeasurementAtGivenTime = round( t/measurementTimeStep) + 1; % + 1 um Index 0 zu vermeiden
p = measuredPressure( indexOfMeasurementAtGivenTime );

function T = temperature(t)

global measurementTimeStep measuredTemperature

% Interpolation möglich aber sinnlos bei so trägen Größen wie Temperatur, Druck und Konzentration
indexOfMeasurementAtGivenTime = round( t/measurementTimeStep) + 1;
T = measuredTemperature( indexOfMeasurementAtGivenTime );

function cox = conOh(t)

global measurementTimeStep measuredconOh
% Interpolation möglich aber sinnlos bei so trägen Größen wie Temperatur, Druck und Konzentration
indexOfMeasurementAtGivenTime = round( t/measurementTimeStep) + 1;
cox = measuredconOh( indexOfMeasurementAtGivenTime );

function coh = conOx(t)

global measurementTimeStep measuredconOx
% Interpolation möglich aber sinnlos bei so trägen Größen wie Temperatur, Druck und Konzentration
indexOfMeasurementAtGivenTime = round( t/measurementTimeStep) + 1;
coh = measuredconOx( indexOfMeasurementAtGivenTime );

Funktion ohne Link?

Gruß

Blues

Chris987

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 12.03.13

	Wohnort: ---

	Version: ---

Verfasst am: 25.04.2013, 21:43 Titel:

Wow! Ich bin schwer beeindruckt.
1. Von der hier gezeigten Kompetenz
2. und noch mehr von der Hilfsbereitsschaft und dem Arbeitseinsatz dafür!

Ein großes Lob und großer Dank! Ich hoffe es der Gemeinschaft zurückgeben zu können. Mal schauen wann ich deinen Level erreiche Smile

Für die Partikelanzahl habe ich in der Tat Messwerte. Zumindest das, was aus rauskommt wird gemessen Wink

Ein Verbrennungsprozess ist nur schwer "live" zu vermessen.

Hab mir auch Tools zur Optimierung gebastelt. Aber noch für die alte Euler-Cauchy-Methode ( der Prof hat es so bezeichnet^^). Das war in der Tat ein Rumspielen an den Konstanten. Deinen Rohling werde ich analysieren und anwenden! Heißen Dank!

Einen gescheiterter Ansatz von mir war wie folgt:
Druck und Temperatur mit interpolationspolynomen verarbeitet. Mit Splines sah es dann auch sehr gut aus. "dsolve" hat sich dann aber geweigert.

Mit besten Grüßen
Christian

PS: wenn ich fragen darf: Was ist dein persönlicher Hintergrund? Informatik?Physik? Mathematik? Oder klassisch Dipl.-Ing? Smile

Bluesmaster

Forum-Century


	Beiträge: 203

	Anmeldedatum: 13.11.11

	Wohnort: Gera

	Version: 2012a

Verfasst am: 25.04.2013, 22:38 Titel:

maschinenbau, deswegen solltest du meiner Lösung nicht zu sehr vertrauen ^^ aber es leitet dich vielleicht in die richtige Richtung.

Im Prinzip muss man aber sagen, dass dein erster Ansatz schon
richtig war. Der ODE-Solver interpoliert lediglich ein bisschen
eleganter mit variabler Schrittweite und so, da kannst
du auch ein wenig an den Toleranzen von "odeset" rumspielen.

Wenn die DGL sehr stark ausschlägt (steif ist) braucht man eventuell
einen anderen Solver. Für dsolve ist die Gleichung eine Nummer
zu groß, da geht es ja um analytische Lösungen.

Das Einzige was mich noch überrascht ist wie man so träge Größen
wie Temperatur und Druck und Konzentration im Mikrosekundenbereich misst?

Gruß

Blues

Chris987

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 12.03.13

	Wohnort: ---

	Version: ---

Verfasst am: 26.04.2013, 09:31 Titel:

Freut mich Herr Kollege Smile

Gemessen werden Luft- und Kraftstoffmassen(und ein bischen mehr). Mit entsprechenden Quarzen kann man den Druck messen. Das ist das entscheidende. Einige andere Werte basieren auf Annahmen (aus der Literatur), die ich an den Summenbrennverlauf gekoppelt habe. Da ich noch einige andere Größen wie das Volumen kenne, kann ich die Gesetze der Thermodynamik anwenden. Thermische Zustandsgleichung in erster Linie. In manchen Fällen sind Druck und Temperatur gar nicht so träge Wink

Über dsolve hat mich "Jörn Loviscach" ein bischen aufgeklärt:) Dachte ich mir dann aber, dass eine analytische Lösung zu viel verlangt ist.

Falls ode45 zu lange braucht, werde ich ode15s testen. Wenn die DGL doch steif sein sollte.

Gruß
Chris


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.