---

Diffusor.tar.gz

Hallo FOAMer,

ich möchte derzeit einen einfachen Diffusor simulieren.
Leider bekomme ich nach einigen Timesteps immer auf eine Gleitkomma Ausnahme:

Zitat:

---

Time = 8.8

smoothSolver:  Solving for Ux, Initial residual = 0.6994309, Final residual = 33.9238507, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.699575075, Final residual = 0.161757026, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.361697662, Final residual = 95.9833167, No Iterations 1000
#0  Foam::error:   rintStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#6  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#7  Foam::fvMatrix<double>::solveSegregated(Foam::Dictionary const&) at ??:?
#8  Foam::fvMatrix<double>::solve(Foam::Dictionary const&) at ??:?
#9  Foam::fvMatrix<double>::solve() at ??:?
#10 
at ??:?
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12 
at ??:?
Floating point exception (core dumped)

---

Und kann nicht wirklich deuten, wo mein Fehler liegt. Ich bin mir noch nicht einmal wirklich sicher wo ich nach dem Fehler suchen soll.
Ich vermute jedoch, dass es eventuell was mit den Anfangsbedingungen im 0-Ordner zu tun haben könnte, denn bei den boundary fields wurde ich aus den Tutorials nicht so richtig schlau. Ich habe zwar versucht, mir das irgendwie aus den Tutorials zusammenzureimen, aber bin mir da als OF-Neuling aber noch recht unsicher.

Desweiteren verstehe ich nicht genau wie ich die turbulent length scale berechnen soll. Dieses ist zB im Programmers Guide als "... If we estimate the turbulent length scale l to be 10% of the width of the inlet... " angenommen. Beziehen sich die 10% dann im 3D-Fall auf den Inletdurchmesser?

Schönen Gruß,

Christoph

[edit: Ich hab euch mal den case mit reingestellt, evtl könnt ihr ihn euch ja mal anschauen... ]

------------------
...im CAD hats gepasst.

[Diese Nachricht wurde von Turbo Encabulator am 09. Dez. 2014 editiert.]

[Diese Nachricht wurde von Turbo Encabulator am 09. Dez. 2014 editiert.]

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Wenn ich den Fehlerbericht richtig deute, verursacht der GaussSeidelSmoother das Problem.
Besonders seltsam finde ich

Zitat:

---

.
.
.

Time = 8

smoothSolver:  Solving for Ux, Initial residual = 0.684048315, Final residual = 0.0442911734, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.997574593, Final residual = 0.0485146775, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 0.699863133, Final residual = 0.0458820356, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.995798716, Final residual = 0.0314446197, No Iterations 1
time step continuity errors : sum local = 5.07312192e+15, global = -4.73147024e+15, cumulative = -4.73134861e+15
smoothSolver:  Solving for epsilon, Initial residual = 0.116120915, Final residual = 4.15924203e-14, No Iterations 1
bounding epsilon, min: -1.02937509e+27 max: 4.62822223e+34 average: 2.38582413e+30
smoothSolver:  Solving for k, Initial residual = 0.478792706, Final residual = 0.0299476502, No Iterations 1
bounding k, min: -6.00531947e+12 max: 6.6653228e+24 average: 2.18093217e+21
ExecutionTime = 179.85 s  ClockTime = 181 s

Time = 8.2

smoothSolver:  Solving for Ux, Initial residual = 0.219633091, Final residual = 0.0123313406, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.306399815, Final residual = 0.0220067635, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 0.218650317, Final residual = 0.0138022874, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.10790366, Final residual = 0.00513502998, No Iterations 1
time step continuity errors : sum local = 2.83931353e+27, global = 1.05148361e+26, cumulative = 1.05148361e+26
smoothSolver:  Solving for epsilon, Initial residual = 0.999999932, Final residual = 0.0239512892, No Iterations 3
bounding epsilon, min: -2.39327517e+31 max: 1.65628566e+55 average: 6.10674604e+51
smoothSolver:  Solving for k, Initial residual = 0.00357014267, Final residual = 4.06804507e-09, No Iterations 1
bounding k, min: -5.1629862e+18 max: 1.95288392e+42 average: 9.48170067e+37
ExecutionTime = 180.09 s  ClockTime = 181 s

Time = 8.4

smoothSolver:  Solving for Ux, Initial residual = 0.590378819, Final residual = 0.0460025547, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.764472835, Final residual = 0.0354139846, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.405278355, Final residual = 0.035982356, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.915739933, Final residual = 0.0855459752, No Iterations 199
time step continuity errors : sum local = 6.68658125e+29, global = 2.60336978e+27, cumulative = 2.70851814e+27
smoothSolver:  Solving for epsilon, Initial residual = 0.0736384258, Final residual = 2.15096757e-21, No Iterations 1
bounding epsilon, min: -2.09754404e+61 max: 4.59707072e+71 average: 4.33773929e+67
smoothSolver:  Solving for k, Initial residual = 0.477796858, Final residual = 0.0300832302, No Iterations 1
bounding k, min: -1.03031232e+38 max: 3.47906758e+58 average: 3.14521036e+54
ExecutionTime = 181.47 s  ClockTime = 182 s

Time 8.6

smoothSolver:  Solving for Ux, Initial residual = 0.00133381417, Final residual = 4.3459851e-05, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.000359973499, Final residual = 3.03707606e-05, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 0.000175138032, Final residual = 9.10573493e-06, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.500223608, Final residual = 0.0085649238, No Iterations 1
time step continuity errors : sum local = 5.4192443e+53, global = -5.24283893e+53, cumulative = -5.24283893e+53
smoothSolver:  Solving for epsilon, Initial residual = 0.15683359, Final residual = 0.00132592472, No Iterations 3
bounding epsilon, min: -2.88228396e+84 max: 6.0033661e+113 average: 1.0096812e+110
smoothSolver:  Solving for k, Initial residual = 2.83687666e-17, Final residual = 2.83687666e-17, No Iterations 0
ExecutionTime = 181.71 s  ClockTime = 183 s

Time 8.8

smoothSolver:  Solving for Ux, Initial residual = 0.6994309, Final residual = 33.9238507, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.699575075, Final residual = 0.161757026, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.361697662, Final residual = 95.9833167, No Iterations 1000
#0  Foam::error:  rintStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?

---

Vorher sind die Initial residuals meist irgendwo im 0.X-Bereich, die Final residuals bei 0.0X.
Dann sackt es im Schritt 8.6 ab, auf Initial 0.000X bzw Final X e-06.
Im Schritt 8.8 "normalisieren" sich die Initials auf 0.X (soweit ich das eben beurteilen kann) aber die Finals heben in X und Z Richtung total ab.
Das gleiche bei No Iterations...

Kann mir jemand auf die Sprünge helfen, warum das so ist?

------------------
...im CAD hats gepasst.

[Diese Nachricht wurde von Turbo Encabulator am 09. Dez. 2014 editiert.]

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Hallo,

deine Deutung des Outputs ist nicht richtig. Du hast ein Problem mit deinem Turbulenzmodell, dass dir dann alles zusammenbrechen lässt. Kennst du pyFoam? Wenn nicht, solltest du dir das mal anschauen. Damit kannst du wunderbar deine Residuen darstellen.

Dein Problem hier ist nun an folgenden Stellen zu suchen:

• Boundary Conditions von U, p, k, epsilon
  
• Numerische Schemen sind ggf. die Ursache
  
• Numerisches Netz (checkMesh)
  
• Schlechte Abschätzung der initialField von den Turbulenzgrößen
  
• Solver settings, in deinem Fall explizit die Unterrelaxation -> SIMPLE
  

Desweiteren berechnet sich die mixtureLength in FOAM, welche in den BC für epsilon und omega verwendet werden sollten wie folgt:

Code:

---

l = d_hyd * 0,0037

---

Der hyraulische Durchmesser lässt sich ja sehr einfach ermitteln.

Anmerkung
Nach der Analyse deines controlDict frage ich mich wieso du dt nicht bei 1 lässt? Ich möchte dir hiermit nur mitteilen, dass dieser Solver keine Zeitableitungen löst und somit das dT keinen Einfluss auf das Ergebnis hat. Die Zeit ist eine Pseudozeit und ist im eigentlichen Sinne die Iterationszahl des SIMPLE Algorithmuses.

Klartext:

dt = 0.005 und du simulierst bis 5s
dt = 0.013 und du simulierst bis 13s
dt = 1 und du simulierst bist 1000s

Die Ergebnisse sind identisch, da die Zeit kein wirkliches Zeitmaß ist. In allen Fällen rechnest du ausschließlich 1000 Iterationen.

Solltest du weitere Hinweise benötigen, kannst du dich ja wieder melden. Viel Erfolg.

------------------
Best regards,
Tobias Holzmann

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Hey danke Tobi!

Hm nur so rein aus Interesse, woran erkenne ich denn jetzt dass der Fehler am Turbulenzmodell liegt?
Wegen dem GAMGSolver im Fehlerbericht?

- Numerisches Netz (checkMesh)
Das Netz sollte eigentlich passen, habe checkMesh durchlaufen lassen, und das sagt es wär ok.

- Solver settings, in deinem Fall explizit die Unterrelaxation -> SIMPLE
Die Unterrelaxation habe ich ein wenig verringert, aber der Fehler bleibt der gleiche. Also denke ich liegt mein Problem nicht daran. Vielleicht kann ich hier später noch ein wenig dran schrauben, wenn der Fehler behoben ist.

Boundary Conditions von U, p, k, epsilon
In meinem Fall habe ich nur Massenstrom, Temperatur und den Druck am Einlass gegeben.
Durch Druck und Temperatur kann ich Dichte des Fluids usw alles ermitteln.
Ich habe mich hier an die Tutorials gehalten, würde mich aber freuen wenn du mal drüberschauen könntest.

- Numerische Schemen sind ggf. die Ursache
Denen werd ich mich jetzt als nächstes widmen...

Schlechte Abschätzung der initialField von den Turbulenzgrößen
Was meinst du damit?

PyFoam funktioniert aber soweit ich weiss nicht mit dem 3.4er Python oder?
Schaut aber auf jeden Fall mal sehr interessant aus, muss ich zugeben... 

------------------
...im CAD hats gepasst.

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Damit auch mal jemand anderes was schreibt, hier mal ein kurzer Hinweis. Soweit ich das weiss, ist es nicht von Vorteil wenn du am Eingang Druck und Massen- bzw. Volumenstrom angibst. Es ist immer besser am Eingang Massenstrom oder Druck anzugeben und am Ausagng enstprechend das andere.

Zum Thema intialFields: Du musst am Anfang deine Trubulenzgrößen abschätzen und angeben. Sind diese zuweit weg von der den "wriklichen" Größen wirst du große Krrekturschritte machen müssen und große Schritte ist immer schlecht da instabil. Ich hoffe meine Angaben stimmen, lasse mich aber gerne berichtigen, so hätte ich auch mal wieder was gelernt.

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Zitat:

---

Original erstellt von User1000:
Damit auch mal jemand anderes was schreibt, hier mal ein kurzer Hinweis. Soweit ich das weiss, ist es nicht von Vorteil wenn du am Eingang Druck und Massen- bzw. Volumenstrom angibst. Es ist immer besser am Eingang Massenstrom oder Druck anzugeben und am Ausagng enstprechend das andere.

---

Wunderbar, das ist wahrhaftig der Fehler... ist mir beim Überfliegen gar nicht aufgefallen.
Du darfst niemals Druck und Geschwindigkeit am Einlass vorgeben sofern dein Solver inkompressible ist. Das liegt an der Natur der Gleichungen. In deinem Fall ist es Überbestimmt und ist damit geneigt in irgendeine Richtung abzudriften. Natürlich kannst du auch den korrekten Wert erzielen. Hängt vom System ab. Vorzustellen ist das wie folgt:

Gauss Elimination
Man ersetzt die untere Dreiecksmatrix mit Nullen, sodass die letzte bekannte Unbekannte bekannt ist. Diese kann man dann für die nächst höhere Unbekannte verwenden usw. Wenn aber die Vorletzte auch schon bekannt ist, ist das System überbestimmt. D.h. mehr Gleichungen als Unbekannte und das stellt oft eine Instabilität dar. Ein begabter Mathematiker würde dir das sicherlich noch wesentlich besser erklären (:

Wie du das siehst das das Turbulenzmodell das Problem ist
Ganz einfach, da deine Turbulenzgrößen utopische Werte annehmen und dann limitiert werden müssen (bounding):

Code:

---

smoothSolver:  Solving for epsilon, Initial residual = 0.116120915, Final residual = 4.15924203e-14, No Iterations 1
bounding epsilon, min: -1.02937509e+27 max: 4.62822223e+34 average: 2.38582413e+30
smoothSolver:  Solving for k, Initial residual = 0.478792706, Final residual = 0.0299476502, No Iterations 1
bounding k, min: -6.00531947e+12 max: 6.6653228e+24 average: 2.18093217e+21

---

Negative Werte sind sowieso unphysikalisch ;) und dann noch -e+12 bzw. -e+24
D.h. diese Werte werden auf Null gesetzt. Das kann hin und wieder beim Beginn der Simulation vorkommen aber sollte nach einigen Iterationen verschwinden.

Zitat:

---

Zum Thema intialFields: Du musst am Anfang deine Trubulenzgrößen abschätzen und angeben. Sind diese zuweit weg von der den "wriklichen" Größen wirst du große Krrekturschritte machen müssen und große Schritte ist immer schlecht da instabil. Ich hoffe meine Angaben stimmen, lasse mich aber gerne berichtigen, so hätte ich auch mal wieder was gelernt.

---

Stimm ich auch zu. Als Anmerkung soll jedoch noch folgendes gegeben werden. Je nach Modell ist die grobe Schätzung mehr oder weniger von Bedeutung. Dem k-Epsilon sagt man ja nicht umsonst hinterher, dass es relativ unempfindlich gegenüber diesen Werten ist. Natürlich sollte das immer in Relation gesehen werden.

PS: Lernst du hier nichts neues  

PPS: Für turbulente Größen am Inlet immer "k » turbulentIntensityKineticEnergyInlet, omega » turbulentMixingLengthDissipationRateInlet, omega » turbulentMixingLengthFrequencyInlet"

Die Namen der BC können auch falsch geschrieben sein (:
Ansonsten eine Abschätzung der Größen kann im hier gemacht werden: Turbulence Estimation

------------------
Best
Tobias Holzmann

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Doch, Doch, lerne ständig was, war ganz anders gemeint, eher nach dem Motto dann hab ich schon wieder was gelernt. Bin jetzt erst aml froh, das meine Aussagen gestimmt haben.

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Selbst wenn se falsch gewesen wäre, hättest dazugelernt (:
Bin auch kein Spezialist und lerne auch immer was dazu... Gibt wesentlich bessere FOAMer - auch hier im Forum - ...

PS: Hab übrigens vor 3 Jahren auch einen Diffusor untersucht 

------------------
Best regards,
Tobias Holzmann

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Hey danke euch zwei!
Ihr seid mir echt eine Große Hilfe!

Habs mal ausprobiert, aber es ist nur eine Iteration weiter gelaufen, bevor der Solver wieder abgehauen ist.
Ich habe jetzt die boundary fields von diff und outlet getauscht.
Also von

Zitat:

---

    outlet
    {
        type            zeroGradient;
    }

    diff
    {
        type            epsilonWallFunction;
        value          uniform 19;
    }

---

auf

Zitat:

---

    outlet
    {
        type            epsilonWallFunction;
        value          uniform 19;
    }

    diff
    {
        type            zeroGradient;
    }

---

Jetzt läuft der case durch. Stimmt das jetzt so, oder hab ich mich da jetzt noch mehr verzettelt?

Tobi, du schreibst l = d_hyd * 0,0037.
Müsste das wenn dann nicht 0,038 heissen? (Also Faktor 10 Unterschied) CFD_Online
Außerdem wäre dem Link nach die Mixting Length 0,07? Woher weiss ich genau, welche von beiden ich hernehmen muss? Welchen benötigt OpenFOAM?
(Habe dazu auf die Schnelle keine weitere Literatur gefunden, die 10% aus dem Tutorials wären jedenfalls näher an der 0,07...)

pyFoam funktioniert bei mir irgendwie noch nicht...

------------------
...im CAD hats gepasst.

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP

---

Hallo Christoph,

bezüglich deiner Randbedingungen für die Geschwindigkeit musst du an den Diffusorwänden auch die Haftbedingung (Ux=Uy=Uz=0) berücksichtigen, eine Null-Gradienten-RB ist nicht korrekt!

Deine Formel l=0.038*dh ist korrekt    

Deine numerische Schemen sind für einen ersten Lauf in Ordnung und so ziemlich das stabilste (1. Ordnung), was OpenFOAM dir bietet.

------------------
Beste Grüße,
Robert

[Diese Nachricht wurde von slint am 11. Dez. 2014 editiert.]

Eine Antwort auf diesen Beitrag verfassen (mit Zitat/Zitat des Beitrags) IP