CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Temperature anomoly at pressure reference cell (http://www.cfd-online.com/Forums/openfoam-solving/84990-temperature-anomoly-pressure-reference-cell.html)

will.logie February 14, 2011 16:17

Temperature anomoly at pressure reference cell
 
1 Attachment(s)
Hello Foamers,

I've been witnessing temperature divergence at the pressure reference cell in a buoyantBoussinesqPimpleFoam problem I've defined axi-symmetrically; a cylindrical hot water tank with a helix-coiled heat exchanger became a 2 wedge of 1 cell thickness. The picture attached hopefully gives an idea of this. In it you can also see where I placed the pressure reference and how the temperature (shown here at 25 seconds) at this point is rising (inexplicably for my liking) and causing all sorts of mischief.

Here my fvSolution:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.7.0                                |
|  \\  /    A nd          | Web:      www.OpenFOAM.com                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p_rgh
    {
        solver                  PCG;
        preconditioner          DIC;
        tolerance              1e-8;
        relTol                  0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol                  0;
    }

    "(U|T|k|epsilon|omega|R)"
    {
        solver                  PBiCG;
        preconditioner          DILU;
        tolerance              1e-6;
        relTol                  0.1;
    }

    "(U|T|k|epsilon|omega|R)Final"
    {
        $U;
        relTol                  0;
    }
}

PIMPLE
{
    momentumPredictor            yes;
    nOuterCorrectors            1;
    nCorrectors                2;
    nNonOrthogonalCorrectors    0;
    //pRefPoint                    ( 0.25 0.13 0 );
    pRefCell                    684;
    pRefValue                  0;
}

relaxationFactors
{
    "(U|T|k|epsilon|omega|R)"                    1;
    "(U|T|k|epsilon|omega|R)Final"                1;
}

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

and fvSchemes:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.7.0                                |
|  \\  /    A nd          | Web:      www.OpenFOAM.com                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default        Euler;
}

gradSchemes
{
    default        Gauss linear;
}

divSchemes
{
    default        none;
    div(phi,U)      Gauss upwind;
    div(phi,T)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega) Gauss upwind;
    div(phi,R)      Gauss upwind;
    div(R)          Gauss linear;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
    default        none;
    laplacian(nuEff,U) Gauss linear uncorrected;
    laplacian((1|A(U)),p_rgh) Gauss linear uncorrected;
    laplacian(kappaEff,T) Gauss linear uncorrected;
    laplacian(DkEff,k) Gauss linear uncorrected;
    laplacian(DepsilonEff,epsilon) Gauss linear uncorrected;
    laplacian(DomegaEff,omega) Gauss linear uncorrected;
    laplacian(DREff,R) Gauss linear uncorrected;
}

interpolationSchemes
{
    default        linear;
}

snGradSchemes
{
    default        uncorrected;
}

fluxRequired
{
    default        no;
    p_rgh          ;
}


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

Any ideas?

Thanks,
Will.

will.logie February 16, 2011 07:31

Anomaly in the sense that...
 
... the temperature at pRefCell oscillates in time from being the lowest (~265K) to the highest (~310K).

In previous cases where I had not yet included the axi-symmetry of the problem (simple 2D extruded mesh) I was not witnessing this in any way. I am thus inclined to think that my error lies therein but am ill equipped to diagnose further. I'd really appreciate some thoughts from you experts out there!

will.logie February 16, 2011 09:33

checkMesh failed?
 
The only thing I can see that seems remotely suspicious is a point in the front wedge slipping outside a threshold (1e-8 meter?). I created the mesh in what I thought was a rather cleanly parameterised gmsh file:

Code:

/*---------------------------------------------------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.7.x                                |
|  \\  /    A nd          | Web:      www.OpenFOAM.com                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
Build  : 1.7.x-da2e560009e6
Exec  : checkMesh
Date  : Feb 16 2011
Time  : 14:25:58
Host  : pc-spf-163
PID    : 14620
Case  : /home/will/OpenFOAM/will-1.7.x/run/ihx/buoyantBoussinesqPimpleFoam/ihx35mmCourseWedge
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:          11948
    internal points:  0
    faces:            23470
    internal faces:  11521
    cells:            5846
    boundary patches: 5
    point zones:      0
    face zones:      0
    cell zones:      1

Overall number of cells of each type:
    hexahedra:    5761
    prisms:        85
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:    0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch              Faces    Points  Surface topology                 
    tesWall            209      418      ok (non-closed singly connected) 
    front              5846    6017    ok (non-closed singly connected) 
    pipeWall            48      96      ok (non-closed singly connected) 
    back                5846    6017    ok (non-closed singly connected) 
    defaultFaces        0        0        ok (empty)                       

Checking geometry...
    Overall domain bounding box (0 0 -0.00567203) (0.324951 0.3 0.00567203)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 1)
    Wedge front with angle 1.00004 degrees
 ***Wedge patch front not planar. Point (0.260332 0.0457038 0.00454412) is not in patch plane by 3.73728e-08 meter.
    Boundary openness (-2.60064e-19 6.7127e-20 6.3524e-19) OK.
    Max cell openness = 3.04681e-16 OK.
    Max aspect ratio = 9.45051 OK.
    Minumum face area = 1.83948e-07. Maximum face area = 0.000137058.  Face area magnitudes OK.
    Min volume = 1.54238e-09. Max volume = 3.80042e-07.  Total volume = 0.000550207.  Cell volumes OK.
    Mesh non-orthogonality Max: 41.1657 average: 8.53724
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.952198 OK.

Failed 1 mesh checks.

End


will.logie February 18, 2011 11:51

... a bigger example
 
1 Attachment(s)
... after 200 seconds.

eugene April 5, 2011 19:52

Well your checkMesh says the front plane of the wedge is non-planar. If you still have the problem try replacing the wedge boundaries with symmetryPlane - if the problem goes away, then the issue is your mesh not being a perfect enough wedge.

will.logie April 6, 2011 04:42

ascii -> binary
 
If one looks inside the average controlDict file there are these two settings which I never really paid attention to before; writeFormat and writePrecision.

Changing writeFormat to binary worked a charm... (thank you Mattijs!)

Checking geometry...
Overall domain bounding box (0 0 -0.00567203) (0.324951 0.5 0.00567203)
Mesh (non-empty, non-wedge) directions (1 1 0)
Mesh (non-empty) directions (1 1 1)
Wedge front with angle 1 degrees
Wedge back with angle 1 degrees
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (1.06916e-19 5.56146e-20 -8.77993e-19) OK.
Max cell openness = 3.05762e-16 OK.
Max aspect ratio = 4.96254 OK.
Minumum face area = 1.22973e-07. Maximum face area = 3.30134e-05. Face area magnitudes OK.
Min volume = 6.95689e-10. Max volume = 7.61277e-08. Total volume = 0.000910615. Cell volumes OK.
Mesh non-orthogonality Max: 42.3953 average: 9.13043
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.529742 OK.

Mesh OK.

xisto April 13, 2011 13:06

After run makeAxialMesh and collapseEdges.

I found a similar problem. however ascii to binary doesn't work for me.

My grid looks fine and the checkMesh before the makeAxialMesh and collapseEdges reports no errors.

Any suggestions?

Thanks

Carlos

Code:



Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:          2900
    internal points:  0
    faces:            5558
    internal faces:  2659
    cells:            1386
    boundary patches: 7
    point zones:      0
    face zones:      0
    cell zones:      0

Overall number of cells of each type:
    hexahedra:    1287
    prisms:        99
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:    0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch              Faces    Points  Surface topology
    front-back          0        0        ok (empty)
    inlet              14      29      ok (non-closed singly connected)
    outlet              14      29      ok (non-closed singly connected)
    simetry            0        0        ok (empty)
    wall                99      200      ok (non-closed singly connected)
    front-back_pos      1386    1500    ok (non-closed singly connected)
    front-back_neg      1386    1500    ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (0 1 -0.360028) (10 3.035 1.36003)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 1)
    Wedge front-back_pos with angle 22.91 degrees
 ***Wedge patch front-back_pos not planar. Point (0.10101 1 0.5) is not in patch plane by 1.28966e-06 meter.
    Boundary openness (5.88638e-18 -1.0902e-15 -5.47392e-16) OK.
    Max cell openness = 2.36954e-16 OK.
    Max aspect ratio = 2.83371 OK.
    Minumum face area = 0.00215678. Maximum face area = 0.241098.  Face area magnitudes OK.
    Min volume = 0.000218165. Max volume = 0.0227451.  Total volume = 8.04848.  Cell volumes OK.
    Mesh non-orthogonality Max: 21.5012 average: 8.88707
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 1.05346 OK.

Failed 1 mesh checks.

End


will.logie April 14, 2011 04:55

Quote:

Originally Posted by xisto (Post 303480)
Code:

Checking geometry...
    Overall domain bounding box (0 1 -0.360028) (10 3.035 1.36003)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 1)
    Wedge front-back_pos with angle 22.91 degrees
 ***Wedge patch front-back_pos not planar. Point (0.10101 1 0.5) is not in patch plane by 1.28966e-06 meter.
    Boundary openness (5.88638e-18 -1.0902e-15 -5.47392e-16) OK.
    Max cell openness = 2.36954e-16 OK.
    Max aspect ratio = 2.83371 OK.
    Minumum face area = 0.00215678. Maximum face area = 0.241098.  Face area magnitudes OK.
    Min volume = 0.000218165. Max volume = 0.0227451.  Total volume = 8.04848.  Cell volumes OK.
    Mesh non-orthogonality Max: 21.5012 average: 8.88707
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 1.05346 OK.

Failed 1 mesh checks.

End


Which axis are you straddling (z)? The wedge has to straddle the axis which you are collapsing equally and it looks to me like you are offset to one side -> Overall domain bounding box (0 1 -0.360028) (10 3.035 1.36003)

23 is a pretty big wedge too. Maybe try <5...?

Which settings did you use for collapseEdges?

xisto April 14, 2011 05:03

Hello Will and thank you for your reply.

I have already tried with 5 degrees and I get the same result. For the collapseEdges I used 1e-8 and 180 degrees

Code:

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:          2900
    internal points:  0
    faces:            5558
    internal faces:  2659
    cells:            1386
    boundary patches: 7
    point zones:      0
    face zones:      0
    cell zones:      0

Overall number of cells of each type:
    hexahedra:    1287
    prisms:        99
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:    0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch              Faces    Points  Surface topology
    front-back          0        0        ok (empty)
    inlet              14      29      ok (non-closed singly connected)
    outlet              14      29      ok (non-closed singly connected)
    simetry            0        0        ok (empty)
    wall                99      200      ok (non-closed singly connected)
    front-back_pos      1386    1500    ok (non-closed singly connected)
    front-back_neg      1386    1500    ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (0 1 0.411235) (10 3.035 0.588765)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 1)
    Wedge front-back_pos with angle 2.49762 degrees
 ***Wedge patch front-back_pos not planar. Point (0.10101 1 0.5) is not in patch plane by 5.03641e-07 meter.
    Boundary openness (6.12585e-19 -3.61958e-17 -8.53818e-16) OK.
    Max cell openness = 2.13747e-16 OK.
    Max aspect ratio = 2.83371 OK.
    Minumum face area = 0.000222607. Maximum face area = 0.0248842.  Face area magnitudes OK.
    Min volume = 2.25174e-05. Max volume = 0.00234757.  Total volume = 0.830702.  Cell volumes OK.
    Mesh non-orthogonality Max: 21.5012 average: 8.88707
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.516369 OK.

Failed 1 mesh checks.


will.logie April 14, 2011 05:30

Hmmm...
 
Hey Carlos,

I've just had a quick look at my makeAxialMesh example and found that this problem you mention persists for me too when I use makeAxialMesh, regardless of binary and writePrecision.

The example where ascii -> binary worked for me (post from April 6th) does not use makeAxialMesh. There I create the wedge in gmsh (translate mesh -1 about symmetryPlane and extrude/revolve 2 about symmetryPlane) and convert that to foam. So I'm not sure I can help you much further with makeAxialMesh - are you using gmsh?

I'm also not sure if you need to straddle the xy-plane equally or not. Maybe someone else can comment on the interpretation from the user manual:

Quote:

For 2 dimensional axi-symmetric cases, e.g. a cylinder, the geometry is specified as a wedge of small angle (e.g. http://www.openfoam.com/docs/user/img/user436x.png) and 1 cell thick running along the plane of symmetry, straddling one of the coordinate planes

xisto April 14, 2011 05:34

I'm using Pointwise for the grid.

Thanks for your help anyway.

What is curious is that the test case that came with makeAxial reports the same error.

jdiorio April 2, 2012 16:39

I'm seeing the same issue, except the mesh isn't axisymmetric. The mesh is a cube of hexahedral cells, looking at a buoyancy flows (solving for piezometric pressure "pd" instead of "p = pd + rhogh"). "pd" is set to "zeroGradient" on all boundaries, and the reference value for pdEqn is supplied by calculating the total pressure:

if (pd.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
pd = p - rhokamb*gh;
}

Everything runs fine when you specify the temperature of the top/bottom boundaries (e.g. "fixedValue") - but if you change one of the boundaries (say to "zeroGradient"), the temperature in the cell that corresponds to "pRefCell" increases, like will.logie showed. I've tried moving the reference cell around, but wherever you move it, that's where you see the temperature anomaly. It's strange because one of the temperature boundaries is still fixed, so a check to something like:

if (T.needReference())
{
Info << "Applying reference temperature " <<
getRefCellValue(T,pRefCell) << " at point " <<
mesh.C()[pRefCell] << endl;
}

doesn't spit out anything, so I can only assume that the TEqn doesn't need a reference cell. Of course, if I set all the temperature boundaries to "zeroGradient" this loop spits out info, and you can see the temperature in that cell getting changed...

How/where could "pRefCell" be entering into the TEqn? Any thoughts?

Thanks.

jdiorio April 5, 2012 16:36

1 Attachment(s)
To try to offer up some more information...

The error comes during the TEqn.solve() step. If you look at the temperature at the cell "pRefCell" you can see it changing wildly after the solution. Changed the top boundary type on T to "fixedGradient" (instead of trying "zeroGradient"), setting the gradient equal to the temperature gradient for the stratification (linear stratification). Bottom boundary still constant value. This appears to fix the problem running in serial, but in parallel the error occurs - the temperature in the cell that corresponds to "pRefCell" becomes unphysical. Attached is a snapshot of what the parallel solution looks like after about 50 iterations.

Any thoughts on how pRefCell could be interfering with the TEqn, even though one of the temperature equation boundaries is a "fixedValue"? This happens in both OF-2.1.0 and OF-1.6-ext (tried both). Tried different solvers for T as well (PCBiCG and GAMG), same deal. Also, when I use a "fixedGradient" at the top, why would it work in serial but not in parallel?

jdiorio April 10, 2012 08:50

Bump. To try and provide some more info, wrote out the location of the pressure reference cell, and then applying:

TEqn.setReference(pRefCell,T.internalField()[pRefCell]);

to see if it's actually setting a temperature reference. Here it would basically set the reference temperature to whatever it was before the TEqn.solve(). Here is an example output - notice how T[pRefCell] is increasing and not equal to what it was before the solve....

Code:

Time = 0.01

Courant Number mean: 0.00783788 max: 0.346089
DILUPBiCG:  Solving for Ux, Initial residual = 0.00193302, Final residual = 4.01574e-10, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.00193302, Final residual = 4.01574e-10, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.00298749, Final residual = 7.85053e-10, No Iterations 3
DICPCG:  Solving for pd, Initial residual = 1, Final residual = 9.99605e-09, No Iterations 313
time step continuity errors : sum local = 2.48078e-09, global = -2.47807e-09, cumulative = -2.47807e-09
DICPCG:  Solving for pd, Initial residual = 0.00911823, Final residual = 9.80415e-09, No Iterations 283
time step continuity errors : sum local = 2.48337e-09, global = -2.47807e-09, cumulative = -4.95615e-09
pRefCell location = (-7.5 -8.9 -19.3)
Setting reference temperature at (-7.5 -8.9 -19.3) to 299.397
T[pRefCell] = 299.397
DILUPBiCG:  Solving for T, Initial residual = 8.09426e-06, Final residual = 8.37418e-09, No Iterations 1
T[pRefCell] = 299.768
ExecutionTime = 50.21 s  ClockTime = 51 s

Time = 0.02

Courant Number mean: 0.000653823 max: 0.393349
DILUPBiCG:  Solving for Ux, Initial residual = 0.952681, Final residual = 1.70106e-10, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 0.952681, Final residual = 1.70104e-10, No Iterations 3
DILUPBiCG:  Solving for Uz, Initial residual = 0.92743, Final residual = 3.8282e-10, No Iterations 3
DICPCG:  Solving for pd, Initial residual = 0.955734, Final residual = 9.13164e-09, No Iterations 307
time step continuity errors : sum local = 1.65392e-09, global = -1.65143e-09, cumulative = -6.60758e-09
DICPCG:  Solving for pd, Initial residual = 0.0156603, Final residual = 9.83983e-09, No Iterations 288
time step continuity errors : sum local = 1.65325e-09, global = -1.65143e-09, cumulative = -8.25901e-09
pRefCell location = (-7.5 -8.9 -19.3)
Setting reference temperature at (-7.5 -8.9 -19.3) to 299.768
T[pRefCell] = 299.768
DILUPBiCG:  Solving for T, Initial residual = 7.12715e-06, Final residual = 2.1977e-09, No Iterations 1
T[pRefCell] = 299.995
ExecutionTime = 85.37 s  ClockTime = 86 s

I've actually tried this specifying some other cell and some other value, which does fix the temperature at that NEW cell, but the value at the old cell still diverges...

jdiorio August 29, 2012 09:39

So I think I found a solution. The strange temperature behavior at pRefCell was not coming from the temperature equation, but from the pressure equation. During the pressure solution, the value of the density perturbation (rhok) at pRefCell would fluctuate. I'm not 100% certain yet, but it's related to the fact that the flux needs to be adjusted AFTER the inclusion of buoyancy. That is, in pEqn.H:

Code:

surfaceScalarField phiU
    (
        (fvc::interpolate(U) & mesh.Sf())

    );

//adjustPhi(phiU, U, pd);

phi = phiU + (rAUf * ((fvc::interpolate(rhok)) * (g & mesh.Sf())));

adjustPhi(phi, U, pd);

The old version of the code I was using did the adjustment prior to the inclusion of the buoyancy term (see comment in code above). Adjusting the flux after the inclusion of buoyancy makes more sense, since the flux that gets adjusted and the flux used in the pressure equation are now the same. Oddly enough it fixes the issue with the density - bit puzzling but this implementation makes more sense to me and fixes the issue I was having. Perhaps the inconsistency in the fluxes causes rhok to be set at pRefCell during the pressure solution?

Yahoo June 5, 2013 18:55

pRefCell on Parallel Simuations
 
Hi
I have developed a solver which works fine on single processor, but diverges when I am running it on multiple processors. For the later, by changing the cell at which pressure is set (pRefCell), my results change. i.e. pRefCell is affecting the results on parallel simulations! Any hints?

Mojtaba.a August 30, 2013 16:58

Dear all,
For those who have the problem about wedge patch which was mentioned above, I want to inform that I solved the problem using Bernhard Gschaider's post in here:

http://www.cfd-online.com/Forums/ope...tml#post387133

  • You have got to change writePrecision value to 12 and rerun the mesh manipulation. after finishing change writePrecision value back to 6.
  • Run checkMesh and problem is gone.
I hope it helps a bit,
Mojtaba

immortality September 1, 2013 11:21

Hi Mojtaba
how writePrecision can help?

Mojtaba.a September 1, 2013 14:29

Quote:

Originally Posted by immortality (Post 449216)
Hi Mojtaba
how writePrecision can help?

Hi Ehsan, :)

Well more values of writePrecision will make OF to write values with more accuracy.
In wedge problems, having a symmetry geometry is one of the main concerns, therefore we are trying to obtain more accuracy in order to have a much better symmetry in our geometry.
Changing writePrecision will help this.

nimasam March 19, 2014 09:11

Quote:

Originally Posted by jdiorio (Post 379311)
So I think I found a solution. The strange temperature behavior at pRefCell was not coming from the temperature equation, but from the pressure equation. During the pressure solution, the value of the density perturbation (rhok) at pRefCell would fluctuate. I'm not 100% certain yet, but it's related to the fact that the flux needs to be adjusted AFTER the inclusion of buoyancy. That is, in pEqn.H:

Code:

surfaceScalarField phiU
    (
        (fvc::interpolate(U) & mesh.Sf())

    );

//adjustPhi(phiU, U, pd);

phi = phiU + (rAUf * ((fvc::interpolate(rhok)) * (g & mesh.Sf())));

adjustPhi(phi, U, pd);

The old version of the code I was using did the adjustment prior to the inclusion of the buoyancy term (see comment in code above). Adjusting the flux after the inclusion of buoyancy makes more sense, since the flux that gets adjusted and the flux used in the pressure equation are now the same. Oddly enough it fixes the issue with the density - bit puzzling but this implementation makes more sense to me and fixes the issue I was having. Perhaps the inconsistency in the fluxes causes rhok to be set at pRefCell during the pressure solution?

another suggestion which works for me was like that:
Code:

fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

          p_rgh = p-rhok*gh;
          p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

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



All times are GMT -4. The time now is 05:13.