[Uyan]

Hi all,

I am new to OpenFOAM software. I ran a simulation of a Jet in a crossflow simulation using simpleFoam. I need to add a passive scalar transportation equation to the simpleFoam solver.

I went through the following example.
http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam

But my problem of including passive scalar transportation equation with a turbulent Schmidt number is slightly different from this example.

So, How do I add a passive scalar transport equation to the simpleFoam solver?

Appreciate any guidance on this issue.

Thank you.

[Uyan]

For anyone who is trying to solve a problem like this i am posting the method i followed.

create a TEqn as of other UEqn and PEqn for the passive scalar.

 {
    volScalarField kappaEff
    (
        "kappaEff",
         (turbulence->nut()/Prt) +(turbulence->nu()/Pr)
     );

    fvScalarMatrix TEqn
    (
      fvm::div(phi, T)
      - fvm::laplacian(kappaEff, T)
    );

    TEqn.relax();     

    TEqn.solve();

}

Then add the T field into createFields.H header file.

    Info<< "Reading field T\n" <<endl;
    volScalarField T
    (
      IOobject
      (
         "T",
         runTime.timeName(),
         mesh,
         IOobject::MUST_READ,
         IOobject::AUTO_WRITE
       ),
      mesh
    );

Then include the TEqn into simpleFoam.C after turbulence correction.

         turbulence->correct();
          #include "TEqn.H"

Then add a line to read turbulent Prandtl number from transportProperties dictionary.

   singlePhaseTransportModel laminarTransport(U, phi);
  

     dimensionedScalar Pr(laminarTransport.lookup("Pr"));

     dimensionedScalar Prt(laminarTransport.lookup("Prt"));

Don't forget to rename simpleFoam.C and the change Make/files to get a new application.

This example was also useful.
https://openfoamwiki.net/index.php/H...ure_to_icoFoam

[Sugajen]

Hi all,

I followed the procedure as posted by Uyan and found that the scalar transport solution diverging. But the scalar is transported throughout the region as expected. What could be the cause of divergence? It does not make sense to me since the solution looks acceptable, it is the divergence bothering me Below is the output from OF, fvSchemes, fvSolution and U BC. Attached are the solver and case files.

All comments and suggestions are welcome

output

Time = 101

smoothSolver:  Solving for Ux, Initial residual = 0.00120042, Final residual = 3.7046e-06, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 0.00120393, Final residual = 3.70971e-06, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.000836939, Final residual = 2.79653e-06, No Iterations 4
GAMG:  Solving for p, Initial residual = 4.31652e-05, Final residual = 8.86733e-07, No Iterations 9
time step continuity errors : sum local = 3.83691e-06, global = 4.17804e-07, cumulative = 0.00012566
GAMG:  Solving for CF, Initial residual = 0.685496, Final residual = 0.513947, No Iterations 1000
ExecutionTime = 211.28 s  ClockTime = 212 s


    From function virtual bool Foam::regIOobject::readIfModified()
    in file db/regIOobject/regIOobjectRead.C at line 291
    
    Re-reading object controlDict from file "/home/ipalab/OpenFOAM/ipalab-4.1/simplefoamexample/system/controlDict"
Time = 102

smoothSolver:  Solving for Ux, Initial residual = 0.00118109, Final residual = 3.65892e-06, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 0.001184, Final residual = 3.66113e-06, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.000818525, Final residual = 2.73043e-06, No Iterations 4
GAMG:  Solving for p, Initial residual = 4.16364e-05, Final residual = 8.8246e-07, No Iterations 9
time step continuity errors : sum local = 3.81498e-06, global = 4.26037e-07, cumulative = 0.000126086
GAMG:  Solving for CF, Initial residual = 0.684755, Final residual = 0.514736, No Iterations 1000
ExecutionTime = 214.31 s  ClockTime = 215 s

fvScheme

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
    div(phi,CF)      bounded Gauss linearUpwind grad(CF);
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
    }

    CF 
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
    };

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        tolerance       1e-08;
        relTol          0.01;
    }

    nuTilda
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        tolerance       1e-08;
        relTol          0.01;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p               1e-5;
        U               1e-5;
        nuTilda         1e-5;
	CF		1e-5;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
	CF		0.3;
    }
    equations
    {
        U               0.7;
        nuTilda         0.7;
    }
}

U

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

internalField   uniform (0.0 0.0 0.1);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value		uniform (0.0 0.0 0.1);
    }

    outlet
    {
        type            zeroGradient;
    }

    xWalls
    {
        type            noSlip;
    }
    
    yWalls
    {
        type            noSlip;
    }

[tas38]

Does your solution diverge using the original application, i.e. without transport of the passive scalar? This is the first thing I would check.

[Uyan]

Sugajen,

I think, how do you calculate diffusion term coefficient DS is the problem here.

Change the DS value and see how it affects the solver performance then you'll see it depends on DS value.

But why do you use a user defined diffusion coefficient in the first place?

If you look at my post #2,you'll see i have used the diffusion coefficient as a function of eddy viscosity.

[Sugajen]

@Uyan I have constant diffusion coefficient. That is why making it user defined. But why would that matter?

@tas38 That's a good point. It does converge with simpleFoam

[Uyan]

Sugen,


In my case i did not relax scalar equation , i simply solved the pressure velocity loop and then using that velocity i solved the scalar.

Even if i had used TEqn.relax() in the above code, i have not under relaxed the field. So it was a redundant statement.

Try without under relaxing. No need to under-relax a passive scalar.