---

Hallo Foamers,

Da bei meiner eigentlichen Rechnung mit rhosimplefoam Probleme mit der Totaltemperatur und darausfolgend mit der Totalenthalpie auftraten, habe ich mir erstmal einen case mit einem einfachen Rohr erstellt, um schneller eine Lösung für mein Problem zu finden. Mein Problem ist, dass meine Totaltemperatur von inlet zu outlet abnimmt, obwohl adiabat(zeroGradient,oder?) gerechnet wird. Leider bin ich immer noch nicht dahinter gekommen was ich falsch mache...Ich hoffe, dass ihr mir da helfen könnt!

Turbulenzmodell: k-omega-SST
Thermotyp: hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas
Hier mal meine 0-files:

Code:

---

p:

dimensions      [1 -1 -2 0 0 0 0];

internalField  uniform 1.0133e5;

boundaryField
{
    OUTLET
    {
        type            zeroGradient;
    }
    WALL
    {
        type            zeroGradient;
    }
    INLET
    {
        type            totalPressure;
        U U;
psi psi;
phi phi;
gamma         1.4;
        p0       uniform 1.0133e5;
value         uniform 1.0133e5;
    }
}

// ************************************************************************* //

U:

dimensions      [0 1 -1 0 0 0 0];

internalField  uniform (0 0 0);

boundaryField
{

    OUTLET
    {
        type            flowRateInletVelocity;
        flowRate        -0.2;     
        value       uniform (0 0 0);
    }
    WALL
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
    INLET
    {
        type            zeroGradient;
    }
}

// ************************************************************************* //

T:

dimensions      [0 0 0 1 0 0 0];

internalField  uniform 298.15;

boundaryField
{

    OUTLET
    {
        type            zeroGradient;
}
    INLET
    {
        type            totalTemperature;
U U;
psi psi;
phi phi;
gamma         1.4;
        T0       uniform 298.15;
value uniform 298.15;
    }
    WALL
    {
        type            zeroGradient;
    }
}

// ************************************************************************* //

k:

dimensions      [0 2 -2 0 0 0 0];

internalField  uniform 1;

boundaryField
{

    OUTLET
    {
        type            zeroGradient;
    }
    INLET
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity      0.05;
        value          uniform 1;
    }
    WALL
    {
        type            compressible::kqRWallFunction;
        value          uniform 1;
    }
}

// ************************************************************************* //

omega:

dimensions      [0 0 -1 0 0 0 0];

internalField  uniform 5000;

boundaryField
{
    OUTLET
    {
        type            zeroGradient;
    }
    INLET
    {
        type            compressible::turbulentMixingLengthDissipationRateInlet;
mixingLength 0.005;
        value          uniform 5000;
    }
    WALL
    {
        type            compressible: megaWallFunction;
        value          uniform 5000;
    }
}

// ************************************************************************* //

alphat und mut sind default.

---

Vielen Dank im Voraus!

Gruß, Eric

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

---

Welche Version von OF verwendest du?

Bei den Versionen vor Of 2.0 bilanziert der Solver rhoSimpleFoam nicht die totale Enthalpie. Der Geschwindigkeitsanteil der Enthalpie bleibt unberücksichtigt.

Im Solver rhoSimpleFoam von OF 2.0 kannst du eine mögliche Implementierung der totalen Enthalpie in hEq.h nachschauen:

{
    fvScalarMatrix hEqn
    (
        fvm::div(phi, h)
      - fvm::Sp(fvc::div(phi), h)
      - fvm::laplacian(turbulence->alphaEff(), h)
    ==
        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
    );

    hEqn.relax();

    hEqn.solve();

    thermo.correct();
}

Viele Grüße

Chrisi

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

---

Hallo.

Danke für deine Antwort. Ich verwende die Windowsversion von OF-2.0! Deswegen versteh ich das ganze ja auch nicht. Vielleicht habe ich an anderer Stelle etwas falsch gemacht? Ich stell nochmal den relevanten Rest rein!

Code:

---

thermophysicalProps:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType      hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>>>>>;

mixture
{
     specie
     {
nMoles 1; //Stoffmenge in mol
molWeight 28.9; //Molmasse in g/mol
     }
     thermodynamics
     {
Cp 1004; //Wärmekapazität cp
Hf 0;        //Schmelzwärme Hf
     }
     transport
     {
As 1.4584e-06; //Sutherlandkoeffizient
Ts 110.33; //Sutherlandtemperatur
     }
}

// ************************************************************************* //

fvSolution:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver           GAMG;
        tolerance        1e-08;
        relTol           0.05;
        smoother GaussSeidel;
        cacheAgglomeration off;
        nCellsInCoarsestLevel 20;
        agglomerator     faceAreaPair;
        mergeLevels      1;
    }

    U
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        nSweeps          2;
        tolerance        1e-7;
        relTol           0.1;
    }

    "(k|omega)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        tolerance       1e-07;
        relTol          0.1;
    }

    h
    {
       solver          PBiCG;
       preconditioner  DILU;
       tolerance       1e-06;
       relTol          0.1;
//solver smoothSolver;
//smoother DILUGaussSeidel;
//nSweeps 2;
//tolerance 1e-10;
//relTol 0.01;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 1;
    rhoMax            rhoMax [1 -3 0 0 0 0 0] 5;
    rhoMin            rhoMin [1 -3 0 0 0 0 0] 0;

    residualControl
    {
        p               1e-9;
        U               1e-6;
        T               1e-7;

        // possibly check turbulence fields
        "(k|omega)" 1e-9;
    }
}

relaxationFactors
{
    p               0.5;
    rho   0.01;
    U               0.7;
    "(k|omega)"     0.9;
    h               0.7;

}

// ************************************************************************* //

fvSchemes:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
grad(U) Gauss linear;
grad(p) Gauss linear;
}

divSchemes
{
    default         none;
  // UEqn
    div(phi,U)          Gauss upwind;
    div((muEff*dev2(T(grad(U)))))   Gauss linear;
div(((muEff*(grad(U)+grad(U).T()))&U)) Gauss linear;
div((((0.666667*muEff)*div(U))*U)) Gauss linear;
    div(U,p)            Gauss upwind;
  // hEqn
    div(phi,h)          Gauss upwind;
  // omegaEqn
    div(phi,omega)    Gauss upwind;
  // epsilonEqn
    div(phi,epsilon)    Gauss upwind;
  // kEqn
div(phi,k)      Gauss upwind;

  // experimental
    div(phid,p)         Gauss upwind;
div(U) Gauss linear;
}

laplacianSchemes
{
    default         none;
  // UEqn
    laplacian(muEff,U)      Gauss linear corrected;
  // hEqn
    laplacian(alphaEff,h)   Gauss linear corrected;
  // pCorrect
    laplacian((rho*rAU),p) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
  // kEqn
    laplacian(DkEff,k) Gauss linear corrected;
  // omegaEqn
    laplacian(DomegaEff,omega) Gauss linear corrected;
  // epsilonEqn
//    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
  // potential
    laplacian(1,p)          Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
div(U,p) upwind phi;

}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p;
}

// ************************************************************************* //

---

Vielleicht entdeckst du ja hier einen Fehler, welchen ich nicht erkenne!  

Vielen Dank im Vorraus!

Gruß, Eric

edit: Bei den divSchemes nicht wundern, da wir noch einen zusätzlichen Term zur Probe in die hEqn eingefügt haben. Hatte aber auch keine Auswirkungen auf das Ergebnis!

[Diese Nachricht wurde von mump0 am 17. Feb. 2012 editiert.]

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

---

Hallo,

Sorry der Solver in OF 2.0 macht doch keine Totale Enthalpiebilanz!:

Diese musst du selbst implementieren.

Eine Möglichkeit ist diese:

{
   
// so muss die H-Gl für die totale Enthalpie aussehen
volScalarField u2 = 0.5 * magSqr(U);
volTensorField gradU = fvc::grad(U);
volTensorField kin1 = turbulence->muEff() * (gradU + gradU.T());
volScalarField divU = fvc::div(U);
volVectorField kin2 = 2.0/3.0 *turbulence->muEff() * divU * U;
volScalarField kin3 = fvc::div(kin2);
volVectorField kin4 = kin1 & U;
volScalarField kin5 = fvc::div(kin4);

fvScalarMatrix hEqn
    (
        fvm::div(phi, h)
      - fvm::Sp(fvc::div(phi), h)
      - fvm::laplacian(turbulence->alphaEff(), h)
     
    ==

      - kin3
      + kin5
      -fvc::div(phi, u2)
    );

    hEqn.relax();
   
    //h.min(200000);

    //h.max(1760000);

    eqnResidual = hEqn.solve().initialResidual();
    maxResidual = max(eqnResidual, maxResidual);

    thermo.correct();

    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
        << max(thermo.T()).value() << endl;

}

Die Konvergenz des Solvers mit dieser hEqn.H ist aber nicht immer super.

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

---

Ok danke dir! Werd es gleich mal implementieren, durchlaufen und dann mal das ergebnis posten.

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

---

beim compilieren sagt er mir folgenden fehler:

hEqn.H: In function 'int main(int, char**)':
hEqn.H:32:5: error: 'eqnResidual' was not declared in this scope
hEqn.H:33:5: error: 'maxResidual' was not declared in this scope

muss quasi noch deklariert werden...aber wie und als was?  

gruß, eric

[Diese Nachricht wurde von mump0 am 21. Feb. 2012 editiert.]

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

---

Ich habe diese Energiegleichung in OF 1.7 implementiert ich denke in OF 2 könnte es so funktionieren:

Zitat:

---

{
   
// so muss die H-Gl für die totale Enthalpie aussehen
volScalarField u2 = 0.5 * magSqr(U);
volTensorField gradU = fvc::grad(U);
volTensorField kin1 = turbulence->muEff() * (gradU + gradU.T());
volScalarField divU = fvc::Div(U);
volVectorField kin2 = 2.0/3.0 *turbulence->muEff() * divU * U;
volScalarField kin3 = fvc::Div(kin2);
volVectorField kin4 = kin1 & U;
volScalarField kin5 = fvc::Div(kin4);

fvScalarMatrix hEqn
    (
        fvm::Div(phi, h)
      - fvm::Sp(fvc::Div(phi), h)
      - fvm::laplacian(turbulence->alphaEff(), h)
     
    ==

      - kin3
      + kin5
      -fvc::Div(phi, u2)
    );

    hEqn.relax();

    hEqn.solve();

    thermo.correct();

}

---

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

---

Hallo,
der solver funktioniert nun. Habe mir mal ein paar Zwischenergebnisse ausgegeben und mit den alten Ergebnissen zu den gleichen Iterationsschritten verglichen. Sieht auf jeden Fall besser aus... Leider ist mir die Rechnung abgeschmiert und nicht schön konvergiert. Hast du schon Erfahrungen was man da in den fvSchemes ändern kann, damit es besser konvergiert?
Werd jetzt erstmal noch eine Rechnung starten mit niedrigeren relaxationsFactors und mal schauen ob es was bringt.

Erstmal noch mal vielen Dank für die hEqn.h ! Wenn es schön konvergiert ist, hoffe ich, dass meine Totaltemperaturdifferenz endlich mal 0 ist 

Gruß, Eric

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

---

In den fvSchemes verwende ich für die zusätzlichen Terme "Gaus linear"

Evtl. hilft es den Unterrelaxationsfaktor für die Enthalpie h etwas herabzusetzen, um die Konvergenz zu verbessern.

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