[Raza Javed]

Hi Alexeym,


First of all, thank you so much for your reply.


Yes, my geometry has three regions:
1. Heaters
2. pipe
3. Solid Box that contains the heaters and Pipe. (For your reference, I am attaching the image of that solid box only). Three slots are the places, where three heaters will be placed, and the hole is Inlet of the pipe, and on the opposite wall of the solid box, there will be outlet.


The mechanism of heat exchange is Conduction.


And the physical processes occur between domains include:


1. Heat transfer from heaters to solid box.
2. Heat transfer from solid box to pipe.


I hope I answered your questions, if NOT, then please let me know.



Thank you

[alexeym]

So, in fact, there are TWO regions: tube with liquid and solid box. Heaters are just areas, where you have energy sources and you add them to the simulation through fvOption. Maybe due to geometry constraints you have decided to split them into separate region.

Could be make your simulation converge without heaters? Start with uniform temperature (293 K), set heaters to (294 K) through boundary conditions on heater-box interfaces. Check if you simulation converges.

Right now you have many points, where your setup can fail, try to reduce them.

[Raza Javed]

Yes, I have tried this in one another geometry, that is slightly different from this one. The only difference is that the pipe is rectangular and straight, without any curves. (pictures are attached for your reference)


In this geometry, three heaters(fvOptions), and infront of them is a straight rectangular pipe. And this simulation converges.


And these all are in the solid box also. (as you can see in the figure)


I tested this geometry without heaters also (without fvOptions), it converges.


I started getting this error when I changed the shape of pipe, from this straight rectangular to curved cylindrical pipe.



Everything else is same for both these geometries.

[Raza Javed]

And I tried to remove the fvOptions and manually assigned the temperature of 294K on the heater-to-box interface,as you suggested, with the original geometry with curved pipe, and got the same error.


My changeDictionaryDict of one of the heaters is attached:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      changeDictionaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

T
{
    internalField   uniform 293;
    defaultFaces
    {
        type            wall;
    }
    boundaryField
    {
        "defaultFaces"
        {
            type            zeroGradient;
        }
        "plate1_to_.*"
        {
            type            compressible::turbulentTemperatureCoupledBaffleMixed;
            Tnbr            T;
            kappaMethod     solidThermo;
            kappaName       none;
            value           uniform 294;
        }
    }
}

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

My controlDict file is attached below:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     chtMultiRegionSimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         20;

deltaT          0.5;

writeControl    timeStep;

writeInterval   1;

purgeWrite      5;

writeFormat     ascii;

writePrecision  7;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

[alexeym]

Since you say, that change of mesh leads to divergence, guess, it is time to see your checkMesh output, fvSchemes, and fvSolution.

Btw, "plate1_to_.*" should be fixedValue, not compressible::turbulentTemperatureCoupledBaffleMix ed, if you would like to fix temperature.

[Raza Javed]

I changed it to fixedValue like this:



ChangeDictionaryDict file

T
{
    internalField   uniform 293;
    defaultFaces
    {
        type            wall;
    }
    boundaryField
    {
        "defaultFaces"
        {
            type            zeroGradient;
        }
        "plate1_to_.*"
        {
            type            fixedValue;
            //Tnbr            T;
            //kappaMethod     solidThermo;
            //kappaName       none;
            value           uniform 294;
        }
    }
}

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

And I started getting the following error:

--> FOAM FATAL ERROR: 
Attempt to cast type fixedValue to type compressible::turbulentTemperatureCoupledBaffleMixed

    From function To& Foam::refCast(From&) [with To = const Foam::compressible::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField; From = const Foam::fvPatchField<double>]
    in file /home/ubuntu/OpenFOAM/OpenFOAM-4.1/src/OpenFOAM/lnInclude/typeInfo.H at line 114.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::compressible::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs() at ??:?
#3  Foam::mixedFvPatchField<double>::evaluate(Foam::UPstream::commsTypes) at ??:?
#4  Foam::mixedEnergyFvPatchScalarField::updateCoeffs() at ??:?
#5  Foam::fvMatrix<double>::fvMatrix(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensionSet const&) at ??:?
#6  Foam::tmp<Foam::fvMatrix<double> > Foam::fv::optionList::operator()<double>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::word const&) at ??:?
#7  ? at ??:?
#8  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9  ? at ??:?
Aborted (core dumped)

And, checkMesh, fvSchemes and fvSolutions are below:


CheckMesh

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

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           31457
    faces:            358864
    internal faces:   355492
    cells:            178589
    faces per cell:   4
    boundary patches: 3
    point zones:      0
    face zones:       0
    cell zones:       5

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

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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                  
    inlet               24       19       ok (non-closed singly connected)  
    outlet              24       19       ok (non-closed singly connected)  
    defaultFaces        3324     1674     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.02 -0.07 -0.059) (-0.01 -0.01 0.011)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (-3.542138e-18 -1.54006e-18 -4.235165e-19) OK.
    Max cell openness = 5.662842e-16 OK.
    Max aspect ratio = 125.2255 OK.
    Minimum face area = 1.112013e-08. Maximum face area = 1.179828e-05.  Face area magnitudes OK.
    Min volume = 1.387383e-12. Max volume = 1.096243e-08.  Total volume = 4.2e-05.  Cell volumes OK.
    Mesh non-orthogonality Max: 88.07634 average: 21.58284
   *Number of severely non-orthogonal (> 70 degrees) faces: 4086.
    Non-orthogonality check OK.
  <<Writing 4086 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 1.815118 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

fvSchemes

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

ddtSchemes
{
    default steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss upwind;
    div(phi,K)      bounded Gauss linear;
    div(phi,h)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss linear;//bounded Gauss upwind;
    div(phi,K)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss linear;//bounded Gauss upwind;
    div(phi,R)      bounded Gauss upwind;
    div(R)          Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default        Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
    p_rgh;
}

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

fvSolution

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    rho
    {
        solver          PCG
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0;
    }

    p_rgh
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;

        smoother         GaussSeidel;

        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator     faceAreaPair;
        mergeLevels      1;

        maxIter          10;
    }

    "(U|h|k|epsilon)"
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }
}

SIMPLE
{
    momentumPredictor on;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       100000;
    rhoMin          rhoMin [1 -3 0 0 0] 0.2;
    rhoMax          rhoMax [1 -3 0 0 0] 2;
}

relaxationFactors
{
    fields
    {
        rho             1;
        p_rgh           0.7;
    }
    equations
    {
        U               0.7;
        h               0.7;
        nuTilda         0.7;
        k               0.7;
        epsilon         0.7;
        omega           0.7;
        "ILambda.*"     0.7;
    }
}

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

[alexeym]

fixedValue should be on both sides not only on plate1 side. Otherwise compressible::turbulentTemperatureCoupledBaffleMix ed from other side tries to cast fisxedValue into compressible::turbulentTemperatureCoupledBaffleMix ed (and fails).

Your fvSchemes are unsuitable for non-orthogonal mesh.

There should be at least several non-orthogonal correctors on your almost 90 degrees non-orthogonal mesh (and let's put aside a question of how you were able to create such low quality mesh for such simple geometry).

[Raza Javed]

Hello Alexeym,


Thank you so much for taking so much time to look deep into my problem.

Originally Posted by alexeym

Your fvSchemes are unsuitable for non-orthogonal mesh.


There should be at least several non-orthogonal correctors on your almost 90 degrees non-orthogonal mesh (and let's put aside a question of how you were able to create such low quality mesh for such simple geometry).

Actually, I am using Salome to create a mesh, and my knowledge in Salome is also very limited, and I simply generated a single mesh for the complete geometry. Here I have one question:


I have seen two approaches,
1. Creating a single mesh for the complete geometry
2. Creating separate mesh for each region (heaters, box, pipe in my case), and then join them to make a final mesh.


But I exactly don't know which one is better. The second option is little complex that's why I tried to keep it simple and used the first option. But, Now, as I know that my mesh quality is very bad, I am thinking to try the second option.


Will that be a good choice? OR There is some other way I should try to improve my mesh?


Thank you

[alexeym]

You were able to create tetrahedral mesh for your domains, so, I think, you can create surface meshes for all regions. You export these surface meshes into STLs, then you can use snappyHexMesh to create hex-dominant meshes for your regions. You can even make mesh orthogonal by switching of snapping and layer addition steps.


Alternatively you can start with orthogonal mesh, use topoSet and STLs of your regions to make cell sets, and then use splitMeshRegions to split mesh into regions.

But you can have converged solution even on rather shitty mesh.

[Raza Javed]

Originally Posted by alexeym

You were able to create tetrahedral mesh for your domains, so, I think, you can create surface meshes for all regions. You export these surface meshes into STLs, then you can use snappyHexMesh to create hex-dominant meshes for your regions. You can even make mesh orthogonal by switching of snapping and layer addition steps.


Alternatively you can start with orthogonal mesh, use topoSet and STLs of your regions to make cell sets, and then use splitMeshRegions to split mesh into regions.

But you can have converged solution even on rather shitty mesh.

Hello Alexeym,


I was able to solve this issue by refining the mesh. Thank you so much for your help.


Now I need to change the fluid flow in my geometry from laminar to Turbulent. But I am having some problem with alphat,


I don't know what is it for? is it necessary to have alphat in our case when modelling the Turbulent flow? because when I include this in my case, it also disturbs other boundary conditions. because I don't know how to use alphat in our simulation? OR what it actually is? and where it is usually used? because I tried to find help relating to this but I couldn't find anything.


I shall be very thankful if you can help me out in this.


Thank you

[alexeym]

alphat is turbulent heat diffusivity, obviously you need it in a simulation with turbulence.

[Raza Javed]

Thank you so much for your reply.


I copied alphat from one of tutorials of Openfoam, and put its content in my system/changeDictionaryDict,where I am putting all of my boundary conditions. My changeDictionaryDict file is given below:

boundary
{
    inlet
    {
        type            patch;
    }
    outlet
    {
        type            patch;
    }
}

T
{
    internalField   uniform 300;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           $internalField;
            
        }

        outlet
        {
            type            zeroGradient;
            value           $internalField;

        }

        "fluid_to_box"
        {
            type            compressible::turbulentTemperatureCoupledBaffleMixed;
            Tnbr            T;
            kappaMethod     fluidThermo;
            value           uniform 300;
        }
    }
}

U
{
    internalField   uniform (0 0 0);

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform (0 0 -1e-4);
        }

        outlet
        {
            type            zeroGradient;
        }
        "fluid_to_box"
        {
            type            noSlip;
        }
    }
}

epsilon
{
    internalField   uniform 0.01;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform 0.01;
        }

        outlet
        {
            type            inletOutlet;
            inletValue      uniform 0.01;
        }

        ".*"
        {
            type            epsilonWallFunction;
            value           uniform 0.01;
        }
    }
}

k
{
    internalField   uniform 0.1;

    boundaryField
    {
        inlet
        {
            type            inletOutlet;
            inletValue      uniform 0.1;
        }

        outlet
        {
            type            zeroGradient;
            value           uniform 0.1;
        }

        ".*"
        {
            type            kqRWallFunction;
            value           uniform 0.1;
        }
    }
}


p_rgh
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            zeroGradient;
            value           uniform 0;
        }

        outlet
        {
            type            fixedValue;
            value           uniform 0;
        }

        ".*"
        {
            type            fixedFluxPressure;
            value           uniform 0;
        }
    }
}

p
{
    internalField   uniform 0;

    boundaryField
    {
        ".*"
        {
            type            calculated;
            value           uniform 0;
        }
    }
}

alphat
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform 0;
        }

        outlet
        {
            type            calculated;
            value           uniform 0;
        }

        "fluid_to_box"
        {
            type            compressible::alphatWallFunction;
            value           uniform 0;
        }


    }
}

nut
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform 0;
        }

        outlet
        {
            type            calculated;
            value           uniform 0;
        }

        "fluid_to_box"
        {
            type            nutkWallFunction;
            value           uniform 0;
        }


    }
}

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

When I RUN the solver, and after 2 3 iterations, it gives me the following error:

--> FOAM FATAL ERROR: 
Maximum number of iterations exceeded

    From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar) const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar) const) const [with Thermo = Foam::hConstThermo<Foam::rhoConst<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy>]
    in file /home/ubuntu/OpenFOAM/OpenFOAM-4.1/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 66.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy> > > >::calculate() at ??:?
#3  Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::rhoConst<Foam::specie> >, Foam::sensibleEnthalpy> > > >::correct() at ??:?
#4  ? at ??:?
#5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6  ? at ??:?
Aborted (core dumped)

In the solver log file, my time step continuity error is very large, may be this error comes because of that.



Log file is attached below:

Time = 0.2


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.9997316, Final residual = 0.08241705, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.9999182, Final residual = 0.07499477, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.999995, Final residual = 0.05218458, No Iterations 2
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.04574217, No Iterations 2
Min/max T:-2937356 1666931
GAMG:  Solving for p_rgh, Initial residual = 0.1649963, Final residual = 0.0009814857, No Iterations 6
time step continuity errors : sum local = 19514.56, global = 1518.633, cumulative = 1518.632
Min/max rho:2 2
DILUPBiCG:  Solving for epsilon, Initial residual = 0.9999722, Final residual = 2.410651e-09, No Iterations 1
DILUPBiCG:  Solving for k, Initial residual = 0.9999915, Final residual = 0.06704835, No Iterations 2

Solving for solid region box
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.03970744, No Iterations 2
Min/max T:-263707.7 90866.11

Solving for solid region plate1
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.02564906, No Iterations 2
Min/max T:299.996 300.0009

Solving for solid region plate2
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.03224829, No Iterations 2
Min/max T:299.9911 300.0009

Solving for solid region plate3
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.02419874, No Iterations 2
Min/max T:299.9997 300.0009

Solving for solid region hot1
DICPCG:  Solving for h, Initial residual = 0.8450844, Final residual = 0.005572691, No Iterations 2
Min/max T:300.0004 300.0031

Solving for solid region hot2
DICPCG:  Solving for h, Initial residual = 0.8459708, Final residual = 0.005387616, No Iterations 2
Min/max T:300.0005 300.003

Solving for solid region hot3
DICPCG:  Solving for h, Initial residual = 0.8490342, Final residual = 0.005782191, No Iterations 2
Min/max T:300.0005 300.003
ExecutionTime = 0.32 s  ClockTime = 0 s

Time = 0.3


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.9484502, Final residual = 0.04230773, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.9291604, Final residual = 0.02407606, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.8690145, Final residual = 0.02865172, No Iterations 2
DILUPBiCG:  Solving for h, Initial residual = 0.9721585, Final residual = 0.01466686, No Iterations 2

I shall be thankful if you can help.


Thank you

[alexeym]

Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.9997316, Final residual = 0.08241705, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.9999182, Final residual = 0.07499477, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.999995, Final residual = 0.05218458, No Iterations 2
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.04574217, No Iterations 2
Min/max T:-2937356 1666931
GAMG:  Solving for p_rgh, Initial residual = 0.1649963, Final residual = 0.0009814857, No Iterations 6
time step continuity errors : sum local = 19514.56, global = 1518.633, cumulative = 1518.632
Min/max rho:2 2
DILUPBiCG:  Solving for epsilon, Initial residual = 0.9999722, Final residual = 2.410651e-09, No Iterations 1
DILUPBiCG:  Solving for k, Initial residual = 0.9999915, Final residual = 0.06704835, No Iterations 2

Solving for solid region box
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.03970744, No Iterations 2
Min/max T:-263707.7 90866.11

Obviously your simulation diverges. Why it diverges, I can't say. Initial/boundary conditions, thermophysical properties, etc etc.

[Raza Javed]

Hi,


If we have incompressible fluid flow, then while writing the boundary conditions for alphat, for the wall boundary type, do we need to mention "incompressible" like below?


type incompressible::alphatJayatillekeWallFunction;
value uniform 0;


Thank you

[alexeym]

It depends on how you introduce incompressibility. In case of CHT solvers, you solve compressible NS with incompressible EOS, so in fact your flow is compressible.

[Raza Javed]

Originally Posted by alexeym

It depends on how you introduce incompressibility. In case of CHT solvers, you solve compressible NS with incompressible EOS, so in fact your flow is compressible.

So, while having the incompressible flow in CHT solvers, we can use compressible wall functions of alphat?

[Raza Javed]

And Can we use other boundary conditions, other than wallFunctions on 'wall' type boundaries?


Thank you

[alexeym]

You should define

while having the incompressible flow in CHT solvers

You can use whatever boundary condition, which correctly describes your process with a given discretisation.

[Raza Javed]

Originally Posted by alexeym

You should define



You can use whatever boundary condition, which correctly describes your process with a given discretisation.

Hi,


Thank you so much for your reply.



I need to ask one thing related to Mesh. I am sorry I am asking this question in this thread, because I didn't get helpful response.



If I have non-orthogonal faces (>70) in my mesh. Can I use that mesh in my simulation? or shall I need to make a new mesh without this error?



Can I correct this problem using nNonOrthoCorrectors in fvSolution?


How many correctors I need to put in fvSolutions if I have 2000 non-orthogonal faces in my mesh?


Thank you.

[alexeym]

Hi,

You can. It would be better to make new mesh with lower non-orthogonality. You should use not only increased value of non-orthogonal correctors, but also specific discretisation schemes.