CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Development of incompressible, temperature dependent solver

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By agustinvo

Reply
 
LinkBack Thread Tools Display Modes
Old   July 6, 2015, 11:07
Default Development of incompressible, temperature dependent solver
  #1
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
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!
agustinvo is offline   Reply With Quote

Old   July 6, 2015, 15:26
Default
  #2
MSF
New Member
 
Join Date: Apr 2014
Posts: 7
Rep Power: 4
MSF is on a distinguished road
Hello,

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

Greetings
MSF is offline   Reply With Quote

Old   July 7, 2015, 02:48
Default
  #3
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
Quote:
Originally Posted by MSF View Post
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 is offline   Reply With Quote

Old   July 7, 2015, 05:40
Default
  #4
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
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?
agustinvo is offline   Reply With Quote

Old   July 7, 2015, 06:46
Default
  #5
Member
 
ali alkebsi
Join Date: Jan 2012
Location: Strasbourg, France
Posts: 82
Rep Power: 6
kebsiali is on a distinguished road
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)
kebsiali is offline   Reply With Quote

Old   July 7, 2015, 07:23
Default
  #6
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
Quote:
Originally Posted by kebsiali View Post
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.
agustinvo is offline   Reply With Quote

Old   July 7, 2015, 07:54
Default
  #7
Member
 
ali alkebsi
Join Date: Jan 2012
Location: Strasbourg, France
Posts: 82
Rep Power: 6
kebsiali is on a distinguished road
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 is offline   Reply With Quote

Old   July 7, 2015, 07:56
Default
  #8
Member
 
ali alkebsi
Join Date: Jan 2012
Location: Strasbourg, France
Posts: 82
Rep Power: 6
kebsiali is on a distinguished road
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
kebsiali is offline   Reply With Quote

Old   July 7, 2015, 08:16
Default
  #9
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
Quote:
Originally Posted by kebsiali View Post
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 View Post
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.
agustinvo is offline   Reply With Quote

Old   July 8, 2015, 04:57
Default
  #10
Member
 
ali alkebsi
Join Date: Jan 2012
Location: Strasbourg, France
Posts: 82
Rep Power: 6
kebsiali is on a distinguished road
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?
kebsiali is offline   Reply With Quote

Old   July 8, 2015, 05:52
Default
  #11
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
Quote:
Originally Posted by kebsiali View Post
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());
agustinvo is offline   Reply With Quote

Old   July 8, 2015, 07:19
Default laplacianFoam with sorce term
  #12
New Member
 
vipin
Join Date: Jun 2015
Posts: 14
Rep Power: 3
vipin1431 is on a distinguished road
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
vipin1431 is offline   Reply With Quote

Old   July 10, 2015, 06:43
Default
  #13
Member
 
ali alkebsi
Join Date: Jan 2012
Location: Strasbourg, France
Posts: 82
Rep Power: 6
kebsiali is on a distinguished road
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
kebsiali is offline   Reply With Quote

Old   July 10, 2015, 09:40
Default
  #14
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
Quote:
Originally Posted by kebsiali View Post
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 is offline   Reply With Quote

Old   July 13, 2015, 10:40
Default
  #15
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
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 is offline   Reply With Quote

Old   July 15, 2015, 08:11
Default
  #16
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
Quote:
Originally Posted by agustinvo View Post
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.
kebsiali likes this.
agustinvo is offline   Reply With Quote

Old   August 6, 2015, 09:08
Default
  #17
Senior Member
 
Agustín Villa
Join Date: Apr 2013
Location: Brussels
Posts: 118
Rep Power: 5
agustinvo is on a distinguished road
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)
agustinvo is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[ANSYS Meshing] Help with element size sandri_92 ANSYS Meshing & Geometry 13 May 26, 2014 09:05
Time dependent solver of unsteady Navier Stokes Mehrez Main CFD Forum 9 May 19, 2014 05:01
divergence detected in amg solver temperature increasing relaxation sweeps Kashif Rashid FLUENT 3 May 3, 2014 06: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


All times are GMT -4. The time now is 10:37.