[makaveli_lcf]

Hello foamers!

I came up with a following problem: during simulation of pouring hot melt into the casting mold I got a temperature values bigger than pouring temperature! My BCs are only fixedValue at the inlet and at the walls of mold (where temperature is lower than at the inlet), and zeroGradient at nozzel boundaries, flat fixed top "freesurface" and at the outlet.

Here are the snapshot, where red color indicates overheated cells. From top comes nozzel, flat top interface is with zeroGradient BC, vertical mold is with constant temperature BC:

Mold_external.jpg
Oveheat1.jpg

So you can see from second picture, that difference in temperature is ~30 K over pouring value, which is not appropriated. The problem arises also from other surfaces with zeroGradient BC. For example from outer wall of submerged part of the nozzle:

nozzelwall.jpg

Where can this problem come from? I don't have any sources of heat, and energy equation contains only energy dissipation parts:

DT/Dt = div(alpha_eff grad(T))

My idea was that error comes from gradient schemes due to non-ortogonality of the mesh. I use following fvSchemes:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     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,R)      Gauss upwind;
    div(R)          Gauss linear;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

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

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

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

Waiting for your suggestions! It is strange, the OF produces error in such simple case.

[alberto]

I'd give a try with limiters on the gradients and laplacian.

The mesh is not orthogonal, but it is not skewed either from the detail you showed.

Please let us know what you find.

Best,

[makaveli_lcf]

Alberto,
I have followed your advice. First please check if I understood you correctly.
So, I modified fvSchemes

gradSchemes
{
    default         Gauss linear;
    grad(T)        limited Gauss linear;
}

laplacianSchemes
{
 .......
// !!    laplacian(kappaEff,T) Gauss linear corrected;
    laplacian(kappaEff,T) Gauss linear limited 0.5;
.......
}

snGradSchemes
{
//  !!    default         corrected;
    default         limited 0.5;
}

And I still get production of high temperature values, see animation I attached.

Over1.gif

It represents same geometry with box slice cutted out to see what is going on
You can see still overheated cell indicated, T error ~10 K.

Waiting for your comments and following advices!

[makaveli_lcf]

I have found that not a diffusion, but convectional part produces wrong values.
So I splitted energy equation following way:

dT1/dt + div(phi, T1) = 0 (1)
dT2/dt = div(alpha_eff * grad(T2)) (2)

Below are the results only transport equation (1) at the left picture, and equations (2), which accounts only for diffusion at the right picture respectively. "dT" value means temperature overestimation value (above pouring temperature). Flow is present in both cases, it's direction is shown with arrows in the picture on the right.

Convection and diffusion splited.jpg

All schemes which I used for such test are without limiting.
So what are your suggestions in that case? There is nothing about limiting divergence nonorthogonality corrections in either User or Programming Guides.
Here I use for equation (1) following divergence scheme:

    div(phi,T1)                  Gauss linearUpwind Gauss linear;

[alberto]

Originally Posted by makaveli_lcf

All schemes which I used for such test are without limiting.
So what are your suggestions in that case? There is nothing about limiting divergence nonorthogonality corrections in either User or Programming Guides.
Here I use for equation (1) following divergence scheme:

    div(phi,T1)                  Gauss linearUpwind Gauss linear;

Hello,

in the first fvSchemes you posted, you were using upwind interpolation, which is limited by definition, since you simply use the upwind cell center value on the face.

If you use linearUpwind, you can additionally limit the gradient as shown here ( http://www.cfd-online.com/Forums/ope...earupwind.html ), using

 div(phi,T1)      Gauss linearUpwind cellLimited  Gauss linear 1;

Best,

[makaveli_lcf]

Alberto!

Thank you very much for such immediate response! I will try your suggestions and report here the results.

[alberto]

Just one additional thought: did you try to re-mesh the part? What mesher did you use to generate the grid that gives you problems?

Best,

[makaveli_lcf]

It was a Gambit, so actually mesh is the third part product, I should not change it...
Btw, I have already tried

div(phi,T1) Gauss linearUpwind cellLimited Gauss linear 1;

and it produces more errors(((

Alberto, what is your opinion regarding reasons, which can cause such problems with a simple scalar transport equation?

Just for additional information, here is output of checkMesh:

Mesh stats
    points:           378191
    faces:            1625937
    internal faces:   1540645
    cells:            642651
    boundary patches: 10
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     227870
    prisms:        138900
    wedges:        0
    pyramids:      1338
    tet wedges:    0
    tetrahedra:    274543
    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
    SEN_Outer           5072     2916     ok (non-closed singly connected)
    Port_Wall           10588    5414     ok (non-closed singly connected)
    SEN_Wall            16692    12887    ok (non-closed singly connected)
    WW_c                22256    22896    ok (non-closed singly connected)
    WW_b                13272    11878    ok (non-closed singly connected)
    WW_a                7734     6059     ok (non-closed singly connected)
    NarrowWall          2782     3024     ok (non-closed singly connected)
    TopWall             2974     2570     ok (non-closed singly connected)
    Outlet              3442     2762     ok (non-closed singly connected)
    Inlet               480      519      ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.686 -1.05 -0.0816556) (0.686 0.803 0.0820117)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-1.55627e-17 2.36049e-16 1.53059e-17) OK.
    Max cell openness = 3.32553e-16 OK.
    Max aspect ratio = 15.485 OK.
    Minumum face area = 1e-06. Maximum face area = 0.000136551.  Face area magnitudes OK.
    Min volume = 8.95047e-10. Max volume = 8.44667e-07.  Total volume = 0.0974925.  Cell volumes OK.
    Mesh non-orthogonality Max: 65.4515 average: 13.9421
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.959158 OK.

Mesh OK.

End

and here's checkMesh -allGeometry -allTopology

Mesh stats
    points:           378191
    faces:            1625937
    internal faces:   1540645
    cells:            642651
    boundary patches: 10
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     227870
    prisms:        138900
    wedges:        0
    pyramids:      1338
    tet wedges:    0
    tetrahedra:    274543
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Topological cell zip-up check OK.
    Face vertices OK.
    Face-face connectivity OK.
  <<Writing 3 cells with with single non-boundary face to set twoInternalFacesCells
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                   Bounding box
    SEN_Outer           5072     2916     ok (non-closed singly connected)   (-0.157456 -0.22 -0.0476851) (0.157458 4.17144e-15 0.0476866)
    Port_Wall           10588    5414     ok (non-closed singly connected)   (-0.153715 -0.22 -0.021562) (0.153715 1.06581e-17 0.021562)
    SEN_Wall            16692    12887    ok (non-closed singly connected)   (-0.0846087 -1.24345e-17 -0.0599263) (0.0846087 0.803 0.0599263)
    WW_c                22256    22896    ok (non-closed singly connected)   (-0.686 -1.05 -0.0380647) (0.686 4.09273e-15 0.0378848)
    WW_b                13272    11878    ok (non-closed singly connected)   (-0.3 -1.05 -0.0766829) (0.3 8.1554e-15 0.0763984)
    WW_a                7734     6059     ok (non-closed singly connected)   (-0.1 -1.05 -0.0816556) (0.1 8.8413e-15 0.0820117)
    NarrowWall          2782     3024     ok (non-closed singly connected)   (-0.686 -1.05 -0.027) (0.686 1.36424e-15 0.027)
    TopWall             2974     2570     ok (non-closed singly connected)   (-0.686 -1.3667e-14 -0.0816556) (0.686 8.8413e-15 0.0820117)
    Outlet              3442     2762     ok (non-closed singly connected)   (-0.686 -1.05 -0.027) (0.686 -1.05 0.027)
    Inlet               480      519      ok (non-closed singly connected)   (-0.0599263 0.803 -0.0599263) (0.0599263 0.803 0.0599263)

Checking geometry...
    Overall domain bounding box (-0.686 -1.05 -0.0816556) (0.686 0.803 0.0820117)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-1.55627e-17 2.36049e-16 1.53059e-17) OK.
    Max cell openness = 3.32553e-16 OK.
    Max aspect ratio = 15.485 OK.
    Minumum face area = 1e-06. Maximum face area = 0.000136551.  Face area magnitudes OK.
    Min volume = 8.95047e-10. Max volume = 8.44667e-07.  Total volume = 0.0974925.  Cell volumes OK.
    Mesh non-orthogonality Max: 65.4515 average: 13.9421
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.959158 OK.
    Min/max edge length = 0.001 0.01938 OK.
    All angles in faces OK.
    Face flatness (1 = flat, 0 = butterfly) : average = 0.999998  min = 0.996924
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 0 average: 2.69124
 ***Cells with small determinant found, number of cells: 3
  <<Writing 3 under-determined cells to set underdeterminedCells

Failed 1 mesh checks.

End

[alberto]

Originally Posted by makaveli_lcf

It was a Gambit, so actually mesh is the third part product, I should not change it...
Btw, I have already tried
and it produces more errors(((

Alberto, what is your opinion regarding reasons, which can cause such problems with a simple scalar transport equation?

Well, just to have a better understanding I would change the mesh, adding boundary layers of hexahedral cells in the problematic zone, so to have a more regular grid there.

Was this mesh working OK in FLUENT if you tried it?

Best,
Alberto

[alberto]

FYI, similar problems on meshes with tets were described her:

http://www.cfd-online.com/Forums/ope...foam-tets.html

Best,

[makaveli_lcf]

Alberto, thank you for reference!

I will comment there to join our effort.

[makaveli_lcf]

A have a following question:

Is there any reason why we solve energy equation in buoyantBoussinesqPisoFoam
after momentum predictor and not after pressure and mass flux correction?

Have a nice weekend!

Best,

Alexander

[alberto]

Originally Posted by makaveli_lcf

Is there any reason why we solve energy equation in buoyantBoussinesqPisoFoam
after momentum predictor and not after pressure and mass flux correction?

I would say it is done to obtain updated thermo data, as a function of h, before solving for the pressure equation.

Best,