Development of incompressible, temperature dependent solver

agustinvo · July 6, 2015, 12:07

Hello,

I want to develope an incompressible solver where the density, heat capacity, thermal conductivity, viscosity... depend on the temperature, so I expect to see the buoyancy effects (Boussinesq is no applicable to my problem, neither the references I have found are clear in relation of its range of applications).

So I want to integrate this in my solver. The properties are already defined as polynomial regressions, I have included the density and the heat capacity inside the derivatives... but since this is an incompressible solver, the pressure equations has to be defined also. I have to include rho as well as in the other equations, but it only gives me huge courant numbers, and finally my simulation crashes (if I run it with Boussinesq, or compressible solvers it does not crash, so it is not the set up but the pressure equation).

Any of you have some experience on this? Thanks in advance!

MSF · July 6, 2015, 16:26

Hello,

did you take into account, that in incompressible flow the equations in OpenFoam are divided through rho?

Greetings

agustinvo · July 7, 2015, 03:48

Quote:

Originally Posted by MSF

Hello,

did you take into account, that in incompressible flow the equations in OpenFoam are divided through rho?

Greetings

Yes. In UEqn and TEqn everything is ok, the problem comes in the pEqn. I am not able to understand all the code, so I don't know where should I include the density.

agustinvo · July 7, 2015, 06:40

Hello,

I think the problem is more related with this:

Code:

Courant Number mean: 0 max: 0
Time = 0

DILUPBiCG:  Solving for Ux, Initial residual = 1.428062e-16, Final residual = 1.428062e-16, No Iterations 0
DILUPBiCG:  Solving for Uy, Initial residual = 2.792809e-18, Final residual = 2.792809e-18, No Iterations 0
DILUPBiCG:  Solving for T, Initial residual = 1, Final residual = 4.894705e-12, No Iterations 1
DICPCG:  Solving for p_rgh, Initial residual = 1.085581e-16, Final residual = 1.085581e-16, No Iterations 0
DICPCG:  Solving for p_rgh, Initial residual = 1.086498e-16, Final residual = 1.086498e-16, No Iterations 0
ExecutionTime = 1.16 s  ClockTime = 2 s

Courant Number mean: 9.354086e-36 max: 2.343958e-32
Time = 0

DILUPBiCG:  Solving for Ux, Initial residual = 5.468248e-20, Final residual = 5.468248e-20, No Iterations 0
DILUPBiCG:  Solving for Uy, Initial residual = 7.420473e-20, Final residual = 7.420473e-20, No Iterations 0
DILUPBiCG:  Solving for T, Initial residual = 0.4145431, Final residual = 2.029354e-12, No Iterations 1
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.009628832, No Iterations 248
DICPCG:  Solving for p_rgh, Initial residual = 2.703377e-05, Final residual = 8.585544e-07, No Iterations 223
ExecutionTime = 2.59 s  ClockTime = 3 s

Courant Number mean: 142726 max: 216211.7
Time = 0

DILUPBiCG:  Solving for Ux, Initial residual = 0.9752576, Final residual = 5.325225e-07, No Iterations 133
DILUPBiCG:  Solving for Uy, Initial residual = 0.755965, Final residual = 8.880502e-07, No Iterations 132
DILUPBiCG:  Solving for T, Initial residual = 0.9999981, Final residual = 5.466545e-08, No Iterations 88
DICPCG:  Solving for p_rgh, Initial residual = 0.9907768, Final residual = 0.009654161, No Iterations 286
DICPCG:  Solving for p_rgh, Initial residual = 0.001860895, Final residual = 9.4345e-07, No Iterations 334
ExecutionTime = 6.48 s  ClockTime = 7 s

Courant Number mean: 94347.71 max: 239519
Time = 0

DILUPBiCG:  Solving for Ux, Initial residual = 0.9781596, Final residual = 1.990685e-07, No Iterations 11
DILUPBiCG:  Solving for Uy, Initial residual = 0.9980529, Final residual = 5.394266e-07, No Iterations 9
DILUPBiCG:  Solving for T, Initial residual = 0.1367406, Final residual = 4.657122e-07, No Iterations 11
DICPCG:  Solving for p_rgh, Initial residual = 0.3061352, Final residual = 0.003004767, No Iterations 309
DICPCG:  Solving for p_rgh, Initial residual = 0.002705098, Final residual = 9.266301e-07, No Iterations 347
ExecutionTime = 8.49 s  ClockTime = 9 s

Courant Number mean: 5.062457e+07 max: 1.529998e+08
Time = 0

DILUPBiCG:  Solving for Ux, Initial residual = 0.3327966, Final residual = 3.049949e-07, No Iterations 17
DILUPBiCG:  Solving for Uy, Initial residual = 0.2931983, Final residual = 9.687052e-07, No Iterations 13
DILUPBiCG:  Solving for T, Initial residual = 0.5618018, Final residual = 8.354536e-07, No Iterations 13
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.009164537, No Iterations 349
DICPCG:  Solving for p_rgh, Initial residual = 0.01243677, Final residual = 9.123202e-07, No Iterations 410
ExecutionTime = 10.82 s  ClockTime = 11 s

while I have specified a deltaTime=0.001. Why is not moving the time?

kebsiali · July 7, 2015, 07:46

give more details
what solver do you use
did you take an existing solver and modify on it
are you sure you are not using a steady state solver (simpleFoam)

agustinvo · July 7, 2015, 08:23

Quote:

Originally Posted by kebsiali

give more details
what solver do you use
did you take an existing solver and modify on it
are you sure you are not using a steady state solver (simpleFoam)

Hello,

I have modified the buoyantBoussinesqPimpleFoam solver, where I include the density rho instead of rhok. Here you can see my equations:

Code:

  // Solve the Momentum equation
    mu=rho*turbulence->nu();
    volScalarField muEff("muEff",mu+rho*turbulence->nut());
    rhoPhi=linearInterpolate(rho*U) & mesh.Sf();
    
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muEff,U)
      - fvc::div(muEff*dev2(fvc::grad(U)().T()))
     ==
        fvOptions(rho, U)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                )*mesh.magSf()
            )
        );

        fvOptions.correct(U);
    }
    }

Code:

volScalarField alphaEff("alphaEff", alpha + turbulenceThermal->alphat());
    //alphaEff=alpha + turbulenceThermal->alphat();
    volScalarField kThermalEff("kThermalEff",alphaEff*Cp*rho);

    rhoCpPhi = linearInterpolate(rho*Cp*U) & mesh.Sf();

    fvScalarMatrix TEqn
    (
        fvm::ddt(rho*Cp,T)
      + fvm::div(rhoCpPhi, T)
      - fvm::laplacian(kThermalEff, T)
     ==
        fvOptions(rho*Cp,T)
    );

    TEqn.relax();

    fvOptions.constrain(TEqn);

    TEqn.solve();

    //radiation->correct();

    fvOptions.correct(T);

Code:

volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();

    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (
            (fvc::interpolate(rho*HbyA) & mesh.Sf())
          + rAUf*fvc::ddtCorr(rho, U, rhoPhi)
        )
      + phig
    );

    fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);

    // Update the fixedFluxPressure BCs to ensure flux consistency
    setSnGrad<fixedFluxPressureFvPatchScalarField>
    (
        p_rgh.boundaryField(),
        (
            phiHbyA.boundaryField()
          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
           *rho.boundaryField()
        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
    );


    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, p_rgh)
        );

        fvOptions.constrain(p_rghEqn);

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            // Calculate the conservative fluxes
            phi = (phiHbyA - p_rghEqn.flux())/fvc::interpolate(rho);

            // Explicitly relax pressure for momentum corrector
            p_rgh.relax();

            // Correct the momentum source with the pressure gradient flux
            // calculated from the relaxed pressure
            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }

    p = p_rgh + rho*gh;

as you can see the equations are ok, and I am using Euler as ddt scheme.

kebsiali · July 7, 2015, 08:54

in UEqn
why dont you use
+ turbulence->divDevRhoReff(U) instead of
- fvm::laplacian(muEff,U)
- fvc::div(muEff*dev2(fvc::grad(U)().T())) ?

in pEqn you should have the following

surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

and

surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
); phiHbyA += phig;

or

surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(HbyA) & mesh.Sf())
+ rAUf*fvc::ddtCorr(rho, U, rhoPhi)
)
+ phig
);

you should have

fvOptions.makeRelative(phiHbyA); not with rho

you should have

setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
)/(mesh.magSf().boundaryField()*rAUf.boundaryField() )
);

you should have

if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - p_rghEqn.flux();

p_rgh.relax();

U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}

kebsiali · July 7, 2015, 08:56

thats just another way of putting the equations and it follows the method used in vof solver

however why the time does not move is a mistry to me

agustinvo · July 7, 2015, 09:16

Quote:

Originally Posted by kebsiali

in UEqn
why dont you use
+ turbulence->divDevRhoReff(U) instead of
- fvm::laplacian(muEff,U)
- fvc::div(muEff*dev2(fvc::grad(U)().T())) ?

I did that because I found no sense in divDevRhoReff in an incompressible solver. Indeed, it is possible to use it.

Quote:

Originally Posted by kebsiali

in pEqn you should have the following

surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

and

surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
); phiHbyA += phig;

or

surfaceScalarField phiHbyA
(
"phiHbyA",
(
(fvc::interpolate(HbyA) & mesh.Sf())
+ rAUf*fvc::ddtCorr(rho, U, rhoPhi)
)
+ phig
);

you should have

fvOptions.makeRelative(phiHbyA); not with rho

you should have

setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField()
- fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
)/(mesh.magSf().boundaryField()*rAUf.boundaryField() )
);

you should have

if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - p_rghEqn.flux();

p_rgh.relax();

U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}

I have tried this, but the problem is still there. I think is related with the no advancing time, since Courant is increasing without control.

kebsiali · July 8, 2015, 05:57

can i ask how and where you defined your variable rho ?
and dont you expect a new source term to appear in the pressure equation due to the varying rho?

agustinvo · July 8, 2015, 06:52

Quote:

Originally Posted by kebsiali

can i ask how and where you defined your variable rho ?
and dont you expect a new source term to appear in the pressure equation due to the varying rho?

rho is defined in the createFields.H file:

Code:

volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
       // mesh
    baseRho+slopeRho*T+slope2Rho*T*T+slope3Rho*T*T*T+slope4Rho*T*T*T*T+slope5Rho*T*T*T*T*T
    );

In relation with the pEqn, I think the equation is ok (I have checked also buoyantPimplefoam and it is quite similar). The term you mean is phig:

Code:

surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

vipin1431 · July 8, 2015, 08:19

hi
I have created the laplacian foam with source term in user directory. how can i use this solver for my mesh? i got problem in compiling this solver plz help..

thanks in advanced

vipin

kebsiali · July 10, 2015, 07:43

the PEqn is based on the fact that the gradient of velocity is equal to zero i.e., continuity
but when the density is a variable; this condition turns into

$\frac {\partial \rho}{\partial t} + \nabla \cdot ( \rho \rm{v} ) = 0,$

actually what i think is:

without this, the solver is within the pimple loop trying to achieve the old condition without marching in time

agustinvo · July 10, 2015, 10:40

Quote:

Originally Posted by kebsiali

the PEqn is based on the fact that the gradient of velocity is equal to zero i.e., continuity
but when the density is a variable; this condition turns into

$\frac {\partial \rho}{\partial t} + \nabla \cdot ( \rho \rm{v} ) = 0,$

actually what i think is:

without this, the solver is within the pimple loop trying to achieve the old condition without marching in time

Hello

thanks for this idea. I will try it now! I will make you know if it is working.

agustinvo · July 13, 2015, 11:40

Hello,

I have included the rhoEqn in my case, and now the time is increasing! But I still have convergence problems, and my Courant number seems to be always equal to zero:

Code:

Starting time loop

Reading surface description:
    z-Plane

Time = 0

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.9973808, Final residual = 7.190542e-23, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.9977049, Final residual = 7.22382e-23, No Iterations 3
DILUPBiCG:  Solving for T, Initial residual = 1, Final residual = 2.407855e-27, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.008707272, Final residual = 8.686427e-08, No Iterations 305
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.9915434, Final residual = 9.755695e-16, No Iterations 540
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 3.53 s  ClockTime = 4 s

Time = 0.0001

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.2423949, Final residual = 4.052451e-16, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.7640606, Final residual = 2.225911e-16, No Iterations 2
DILUPBiCG:  Solving for T, Initial residual = 0.4897039, Final residual = 1.178537e-27, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 9.357006e-06, No Iterations 365
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.6779252, Final residual = 8.230382e-15, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 7.06 s  ClockTime = 7 s

Time = 0.0002

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.4997963, Final residual = 2.599016e-16, No Iterations 72
DILUPBiCG:  Solving for Uy, Initial residual = 0.8234894, Final residual = 4.587841e-16, No Iterations 81
DILUPBiCG:  Solving for T, Initial residual = 0.9999986, Final residual = 5.121843e-16, No Iterations 16
DICPCG:  Solving for p_rgh, Initial residual = 0.9987744, Final residual = 9.690122e-06, No Iterations 418
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.4040092, Final residual = 0.2838411, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 11.66 s  ClockTime = 12 s

Time = 0.0003

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.7760858, Final residual = 2.422542e-16, No Iterations 33
DILUPBiCG:  Solving for Uy, Initial residual = 0.729375, Final residual = 7.134555e-16, No Iterations 33
DILUPBiCG:  Solving for T, Initial residual = 0.4481602, Final residual = 2.430657e-16, No Iterations 36
DICPCG:  Solving for p_rgh, Initial residual = 0.9999992, Final residual = 9.311596e-06, No Iterations 738
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.5479091, Final residual = 31.94313, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 16.75 s  ClockTime = 17 s

Time = 0.0004

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.9067029, Final residual = 4.97392e-16, No Iterations 9
DILUPBiCG:  Solving for Uy, Initial residual = 0.7489127, Final residual = 3.218888e-18, No Iterations 10
DILUPBiCG:  Solving for T, Initial residual = 0.5418659, Final residual = 1.234434e-16, No Iterations 15
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 1810.988, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.7886131, Final residual = 286.9007, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 22.1 s  ClockTime = 22 s

Time = 0.0005

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.9805554, Final residual = 8.929525e-17, No Iterations 9
DILUPBiCG:  Solving for Uy, Initial residual = 0.9615684, Final residual = 3.209775e-16, No Iterations 8
DILUPBiCG:  Solving for T, Initial residual = 0.9894829, Final residual = 7.77986e-16, No Iterations 10
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 1199.078, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 7.390336e-05, Final residual = 1.163134, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 27.23 s  ClockTime = 27 s

Time = 0.0006

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.6522411, Final residual = 1.133062e-18, No Iterations 8
DILUPBiCG:  Solving for Uy, Initial residual = 0.8703662, Final residual = 7.707482e-17, No Iterations 7
DILUPBiCG:  Solving for T, Initial residual = 0.1164554, Final residual = 1.04748e-16, No Iterations 10
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 8.014764, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 1.679593e-05, Final residual = 494607.5, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 32.36 s  ClockTime = 33 s

Time = 0.0007

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  
 at sigaction.c:?
#3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  
 at ??:?
#6  
 at ??:?
#7  
 at ??:?
#8  
 at ??:?
#9  __libc_start_main at ??:?
#10  
 at ??:?
Floating point exception

There is still something not working properly over there. It seems tobe the pEqn for me:

Code:

{
    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rAUf("rAUf", fvc::interpolate(rho*rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();

    surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rho)*mesh.magSf());

    surfaceScalarField phiHbyA
    (
        "phiHbyA",
        (
            (fvc::interpolate(rho*HbyA) & mesh.Sf())
          + rAUf*fvc::ddtCorr(rho, U, rhoPhi)
        )
      + phig
    );

    fvOptions.makeRelative(phiHbyA);

    // Update the fixedFluxPressure BCs to ensure flux consistency
    setSnGrad<fixedFluxPressureFvPatchScalarField>
    (
        p_rgh.boundaryField(),
        (
            phiHbyA.boundaryField()
          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
           *rho.boundaryField()
        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
    );

//    Info<<"phiHbyA" << phiHbyA<< endl;
//    
//    Info<<"rAUf" << rAUf<< endl;
//    Info<<"p_rgh" << p_rgh<< endl;

    fvScalarMatrix p_rghDDtEqn
    (
        //fvc::ddt(rho) //+ psi*correction(fvm::ddt(p_rgh))
        fvc::div(phiHbyA)
      - fvm::laplacian(rAUf, p_rgh)
 //     ==
 //       fvOptions(p_rgh)
    );

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            p_rghDDtEqn
          
        );

        fvOptions.constrain(p_rghEqn);

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            // Calculate the conservative fluxes
            rhoPhi = phiHbyA - p_rghEqn.flux();

            // Explicitly relax pressure for momentum corrector
            p_rgh.relax();

            // Correct the momentum source with the pressure gradient flux
            // calculated from the relaxed pressure
            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }
    
    p = p_rgh + rho*gh;



    #include "densityEqn.H"
    #include "continuityErrs.H"
}

I will take a look to the setup of the case, but if previously it worked (in buoyantBoussinesqPimpleFoam), I don't know where the problem could be.

agustinvo · July 15, 2015, 09:11

Quote:

Originally Posted by agustinvo

Hello,

Code:

Starting time loop

Reading surface description:
    z-Plane

Time = 0

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.9973808, Final residual = 7.190542e-23, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.9977049, Final residual = 7.22382e-23, No Iterations 3
DILUPBiCG:  Solving for T, Initial residual = 1, Final residual = 2.407855e-27, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 0.008707272, Final residual = 8.686427e-08, No Iterations 305
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.9915434, Final residual = 9.755695e-16, No Iterations 540
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 3.53 s  ClockTime = 4 s

Time = 0.0001

Courant Number mean: 0 max: 0
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Ux, Initial residual = 0.2423949, Final residual = 4.052451e-16, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.7640606, Final residual = 2.225911e-16, No Iterations 2
DILUPBiCG:  Solving for T, Initial residual = 0.4897039, Final residual = 1.178537e-27, No Iterations 2
DICPCG:  Solving for p_rgh, Initial residual = 1, Final residual = 9.357006e-06, No Iterations 365
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
DICPCG:  Solving for p_rgh, Initial residual = 0.6779252, Final residual = 8.230382e-15, No Iterations 1001
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
ExecutionTime = 7.06 s  ClockTime = 7 s

I have found my error: since now my flux is defined as rhoPhi, phi was never updated, so it was always zero, and I was not able to adapt timeStep, neither estimate the Courant number. No I can see it without problem, and I have reduce my timeStep.

It only remains to see if it is going ok now.

In relation with the rhoEqn.H, I have disabled it.

agustinvo · August 6, 2015, 10:08

Hello,

first of all, I was able to create the solver!

Now, I am dealing with the set up of the case. Base in compressible equations (I am simulating air), I should expect some turbulence in my case. The problem comes when, by using the same BC for T and U, there is not!

If I use Boussinesq, I am able to see the turbulence, but temperature differences are long (we speak about 100K). In the other case, by using temperature dependence, the differences are around 20K.

Does anyone up to when Boussinesq is a good approximation? (not for my case, but just to know, since there are not general indications)

July 8, 2015, 08:19	laplacianFoam with sorce term	#12
vipin1431 New Member vipin Join Date: Jun 2015 Posts: 14 Rep Power: 11	hi I have created the laplacian foam with source term in user directory. how can i use this solver for my mesh? i got problem in compiling this solver plz help.. thanks in advanced vipin

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Time dependent solver of unsteady Navier Stokes	Mehrez	COMSOL	10	March 20, 2024 04:57
[ANSYS Meshing] Help with element size	sandri_92	ANSYS Meshing & Geometry	14	November 14, 2018 08:54
divergence detected in amg solver temperature increasing relaxation sweeps	Kashif Rashid	FLUENT	3	May 3, 2014 07:23
icoUncoupledKinematicParcelFoam: temperature dependent viscosity	90nash	OpenFOAM Running, Solving & CFD	5	March 27, 2014 00:43
Temperature dependent propertys and finalIteration in chtMultiRegionFoam	waiter120	OpenFOAM	0	February 20, 2013 06:22

July 6, 2015, 12:07	Development of incompressible, temperature dependent solver	#1
agustinvo Senior Member Agustín Villa Join Date: Apr 2013 Location: Alcorcón Posts: 314 Rep Power: 15	Hello, I want to develope an incompressible solver where the density, heat capacity, thermal conductivity, viscosity... depend on the temperature, so I expect to see the buoyancy effects (Boussinesq is no applicable to my problem, neither the references I have found are clear in relation of its range of applications). So I want to integrate this in my solver. The properties are already defined as polynomial regressions, I have included the density and the heat capacity inside the derivatives... but since this is an incompressible solver, the pressure equations has to be defined also. I have to include rho as well as in the other equations, but it only gives me huge courant numbers, and finally my simulation crashes (if I run it with Boussinesq, or compressible solvers it does not crash, so it is not the set up but the pressure equation). Any of you have some experience on this? Thanks in advance!

July 6, 2015, 16:26		#2
MSF New Member Join Date: Apr 2014 Location: Germany Posts: 24 Rep Power: 12	Hello, did you take into account, that in incompressible flow the equations in OpenFoam are divided through rho? Greetings

July 7, 2015, 07:46		#5
kebsiali Member ali alkebsi Join Date: Jan 2012 Location: Strasbourg, France Posts: 82 Rep Power: 14	give more details what solver do you use did you take an existing solver and modify on it are you sure you are not using a steady state solver (simpleFoam)

July 7, 2015, 08:54		#7
kebsiali Member ali alkebsi Join Date: Jan 2012 Location: Strasbourg, France Posts: 82 Rep Power: 14	in UEqn why dont you use + turbulence->divDevRhoReff(U) instead of - fvm::laplacian(muEff,U) - fvc::div(muEffdev2(fvc::grad(U)().T())) ? in pEqn you should have the following surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); and surfaceScalarField phiHbyA ( "phiHbyA", (fvc::interpolate(HbyA) & mesh.Sf()) + fvc::interpolate(rhorAU)fvc::ddtCorr(U, phi) ); phiHbyA += phig; or surfaceScalarField phiHbyA ( "phiHbyA", ( (fvc::interpolate(HbyA) & mesh.Sf()) + rAUffvc::ddtCorr(rho, U, rhoPhi) ) + phig ); you should have fvOptions.makeRelative(phiHbyA); not with rho you should have setSnGrad<fixedFluxPressureFvPatchScalarField> ( p_rgh.boundaryField(), ( phiHbyA.boundaryField() - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField()) )/(mesh.magSf().boundaryField()rAUf.boundaryField() ) ); you should have if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA - p_rghEqn.flux(); p_rgh.relax(); U = HbyA + rAUfvc::reconstruct((phig - p_rghEqn.flux())/rAUf); U.correctBoundaryConditions(); fvOptions.correct(U); }

July 7, 2015, 08:56		#8
kebsiali Member ali alkebsi Join Date: Jan 2012 Location: Strasbourg, France Posts: 82 Rep Power: 14	thats just another way of putting the equations and it follows the method used in vof solver however why the time does not move is a mistry to me

July 8, 2015, 05:57		#10
kebsiali Member ali alkebsi Join Date: Jan 2012 Location: Strasbourg, France Posts: 82 Rep Power: 14	can i ask how and where you defined your variable rho ? and dont you expect a new source term to appear in the pressure equation due to the varying rho?

July 10, 2015, 07:43		#13
kebsiali Member ali alkebsi Join Date: Jan 2012 Location: Strasbourg, France Posts: 82 Rep Power: 14	the PEqn is based on the fact that the gradient of velocity is equal to zero i.e., continuity but when the density is a variable; this condition turns into $\frac {\partial \rho}{\partial t} + \nabla \cdot ( \rho \rm{v} ) = 0,$ actually what i think is: without this, the solver is within the pimple loop trying to achieve the old condition without marching in time

August 6, 2015, 10:08		#17
agustinvo Senior Member Agustín Villa Join Date: Apr 2013 Location: Alcorcón Posts: 314 Rep Power: 15	Hello, first of all, I was able to create the solver! Now, I am dealing with the set up of the case. Base in compressible equations (I am simulating air), I should expect some turbulence in my case. The problem comes when, by using the same BC for T and U, there is not! If I use Boussinesq, I am able to see the turbulence, but temperature differences are long (we speak about 100K). In the other case, by using temperature dependence, the differences are around 20K. Does anyone up to when Boussinesq is a good approximation? (not for my case, but just to know, since there are not general indications)