[Raza Javed]

Hello Everyone,


I am using chtMultiRegionSimpleFoam and my Openfoam version is 4.1.


I made my geometry in Salome and I am modelling a laminar fluid flow.


Now I want to change the flow to Turbulent, But I don't know what changes do I need to make in my case to model a Turbulent flow.


Can someone please guide me how to do this?


Thank you

[mrishi]

Hi. I am myself new to the platform but I can give you a few pointers.



Start with modifying the constant/turbulenceProperties where you set simulation type from laminar to RAS (or LES). Then there are a few keywords associated with the type of simulation you choose. eg. for RAS you have to specify the turbulence model used, eg kEpsilon.... you can easily search online for detailed list of available models and what are the corresponding keywords.


Now, assuming you use kEpsilon model, the k and epsilon fields will be solved in order to obtain correct effective viscosity in our momentum equation (UEqn). So in your 0 directory, setup two files, named k and epsilon and define initial and boundary conditions. (You may have to declare k and omega if you use SSTkOmega model) These are all scalar quantities so you can start with the pressure file and modify it. Keep in mind to define dimensions etc properly. Alternatively, you can take a look at the high reynolds number cavity flow tutorial where you will find a sample of these files which you can change.



Next, you will have to specify the following fluxes in your fvSchemes dictionary: div(phi,k), div(phi,epsilon) for which you can use Gauss Linear or Gauss vanLeer schemes as a starting point. This will tell the solver how these two scalars interact with the momentum flux in the discretized form of equations.



Finally, you will need to setup numerical solution for the two discretized equations and error tolerance criteria. This is in the fvSolution dictionary where you may use the prescribed method for pressure and use it for the variables k and epsilon. You will be able to find information on this. Try to run this setup now. OF will tell you if it expects additional things to be declared, along with the dictionary where they are expected, which will give you an idea of where to make the changes. For example, you may need to define things like epsilonFinal which is the variable epsilon but used in the final iteration of a timestep (read about the PIMPLE iteration procedure to get an idea about this aspect of inner and outer iterations).



Hope this helps.

Note: at the time of running OF might ask you to also specify a third variable called nut which is not used by the simulation but seems necessary for the framework to work (I may be wrong here, this is just what I experienced in my case).

[Raza Javed]

Thank you so much for your reply.


I tell you what I have done.


First, as you said, I changed my turbulenceProperties file from laminar to the following:


turbulenceProperties

simulationType RAS;


RAS
{
    RASModel            kEpsilon;

    turbulence          on;
    printCoeffs         on;
}

After putting this, I changed my k and epsilon files. But, one thing here is that I have changeDictionaryDict file in each region. So I went to system/fluid/changeDictionaryDict to change the k and epsilon entries.


And my changeDictionaryDict is given below:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      changeDictionaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
boundary
{
    inlet
    {
        type            patch;
    }
    outlet
    {
        type            patch;
    }

}
U
{
    internalField   uniform (0 1e-3 0);

    boundaryField
    {
        inlet
        {
            type                fixedValue;//flowRateInletVelocity;//pressureInletVelocity;//
            //volumetricFlowRate  0.2;
            //extrapolateProfile  yes;
            value               uniform (0 2 0);
        }

        outlet
        {
            type                inletOutlet;//zeroGradient;////
            value               $internalField;//uniform (0 0 0);//
            inletValue          $internalField;
        }
        "fluid_to_box"
        {
            type                noSlip;
        }
    }
}

T
{
    internalField   uniform 300;

    boundaryField
    {
        inlet
        {
            type            fixedValue;
            value           uniform 300;//$internalField;
            
        }

        outlet
        {
            type            inletOutlet;
            value           $internalField;
            inletValue      $internalField;
        }

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


epsilon
{
    internalField   uniform 0.01;

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

        outlet
        {
            type            zeroGradient;//inletOutlet;
            //inletValue      uniform 0.01;
        }

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

k
{
    internalField   uniform 0.1;

    boundaryField
    {
        inlet
        {
            type            fixedValue;//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            fixedFluxPressure;//zeroGradient;
            value           uniform 0;
        }

        outlet
        {
            type            fixedValue;
            value           uniform 0;
        }

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

/*p
{
    internalField   uniform 0;

    boundaryField
    {
        ".*"
        {
            type            zeroGradient;
            //value           uniform 0;
        }
    }
}*/

p
{
    internalField   uniform 0;

    boundaryField
    {
        inlet
        {
            type            zeroGradient;
            //value           uniform 1;
        }
        outlet  
        {
            type            fixedValue;
            value           uniform 0;
        }        
        "fluid_to_.*"
        {
            type            zeroGradient;
        }
    }
}

Then I changed the fvschemes dictionary, and that is also attached below:

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;
}

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

I couldn't understand your last step, relating to fvSolutions. My fvSolutions are given below:

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;
    }
}

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

After changing up to this step, I am getting the following error.

Time = 0.1


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 0.06324153, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.01731178, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 1, Final residual = 0.05947683, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.9998487, Final residual = 0.01933637, No Iterations 2
Min/max T:300 300
GAMG:  Solving for p_rgh, Initial residual = 1, Final residual = 0.7186809, No Iterations 10
time step continuity errors : sum local = 0.8787344, global = -0.03886439, cumulative = -0.03886439
Min/max rho:2 2
DILUPBiCG:  Solving for epsilon, Initial residual = 0.2135582, Final residual = 0.01602102, No Iterations 1
bounding epsilon, min: -66.96574 max: 1520.873 average: 78.71354
DILUPBiCG:  Solving for k, Initial residual = 1, Final residual = 0.02653225, No Iterations 2

Solving for solid region box
DICPCG:  Solving for h, Initial residual = 0.8551727, Final residual = 0.03800197, No Iterations 2
Min/max T:300 300

Solving for solid region plate1
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.009459444, No Iterations 2
Min/max T:300 300.0969

Solving for solid region plate2
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.00875892, No Iterations 2
Min/max T:300 300.0972

Solving for solid region plate3
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.008858153, No Iterations 2
Min/max T:300 300.1049
ExecutionTime = 1.56 s  ClockTime = 1 s

Time = 0.2


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.6605907, Final residual = 0.02840858, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.6064321, Final residual = 0.0232098, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.5425838, Final residual = 0.0368132, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.1390125, Final residual = 0.003045702, No Iterations 2
Min/max T:-3.050845e+07 3.376863e+07
GAMG:  Solving for p_rgh, Initial residual = 0.1220398, Final residual = 0.0010412, No Iterations 7
time step continuity errors : sum local = 737166.4, global = 3599.479, cumulative = 3599.44
Min/max rho:2 2
DILUPBiCG:  Solving for epsilon, Initial residual = 0.0009370397, Final residual = 3.433172e-05, No Iterations 1
bounding epsilon, min: -3.263126e+09 max: 1.011055e+10 average: 3.31966e+07
DILUPBiCG:  Solving for k, Initial residual = 0.426422, Final residual = 0.03064464, No Iterations 1
bounding k, min: -4.651673e+07 max: 1.451633e+08 average: 968806.9

Solving for solid region box
DICPCG:  Solving for h, Initial residual = 1, Final residual = 0.03428405, No Iterations 2
Min/max T:-3173387 3195540

Solving for solid region plate1
DICPCG:  Solving for h, Initial residual = 0.9985548, Final residual = 0.01910019, No Iterations 2
Min/max T:-259.3801 682.9494

Solving for solid region plate2
DICPCG:  Solving for h, Initial residual = 0.9995806, Final residual = 0.02288655, No Iterations 2
Min/max T:-7493.021 3313.353

Solving for solid region plate3
DICPCG:  Solving for h, Initial residual = 0.9994745, Final residual = 0.02052039, No Iterations 2
Min/max T:-4433.842 12260.43
ExecutionTime = 2.37 s  ClockTime = 2 s

Time = 0.3


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.9826942, Final residual = 0.0005316068, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.09963916, Final residual = 0.00187412, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.3857384, Final residual = 0.003870226, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.02362706, No Iterations 2


--> 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)

I tried to explain my case as clear as I can. if you need any further clarification please let me know.


Thank you

[mrishi]

I saw your query about your error on multiple posts. I am also curiously following the discussion.
In principle,the turbulence case setup correctly. However, in addition to the consistent error of max iterations exceeded, we can see that the value of k, epsilon are blowing up. In my case I had tackled this by setting boundary conditions correctly and it seems that the BCs you have set up are very loosely defined as almost all of them begin to get unbounded almost immediately.

For starters:
1. Try to figure out the boundary conditions for k and epsilon at inlet by calculating approximate Reynolds number, turbulent mixing length, turbulence intesity ratio (There are prescriptions you can find for simple geometries which you can use.)
2. I am not too sure about the inletOutlet condition for velocity. Try a simpler definition like zeroGradient (atleast for the purpose of stability solution. You could try out better options later). In general you could replace a lot of BCs with simpler conditions like fixedValue of Gradient as these are primitive types of BCs in CFD. This is related to the advice you received about the pressure-velocity boundary conditions. It might also explain high value of residuals as system is unphysical.

3. I expect that with carefully configuring boundaries you will be able to make your solution stable, it takes some time but it is worthwhile to read about standard and suitable types for each variable. Now, coming to the final error, which seems to be related to species eqn, this post might be helpful :
FOAM FATAL ERROR Maximum number of iterations exceeded
^but #3 can be expected to be heavily related to the other points as the divergence appears to propagate across variables before your solution crashes.

These are only some of my hunches. Do let me know if you figure out what's happening.

[Raza Javed]

Thank you so much for your reply.


Here, I would like to add something.



1. Before I had one simple geometry(picture attached, narrow_pipe10), in which there is a rectangular fluid region, and the fluid flow was laminar. with the same boundary conditions. The solver was running upto steady-state.


2. Then, I just changed the shape of fluid region from rectangular to cylindrical, and put some curves in it (picture attached, inside_geometry).


The three plates in both geometries are three heat sources which I made using fvOptions.


My questions here are the following:


1. Does Laminar flow cannot be modeled in this curve shape fluid region? Should I change the flow from laminar to turbulent for this new geometry? because in this new geometry I am getting this error "maximum number of iterations exceeded", with the same boundary conditions as I had for the old geometry.


2. Can we model turbulent flow in chtMultiRegionSimpleFoam?, because I am using this solver.


3. How can I figure out the boundary conditions for k and epsilon using Reynolds number? (I am sorry for this question, as I am very new to Openfoam so I don't know much about it)



These are my first basic confusion here, with same boundary condition one geometry is working while other is NOT. So, I thought may be I need to change the flow from Laminar to Turbulent for this new geometry.


I shall be very thankful if you can clear these doubts.


In the meanwhile, I look for the boundary conditions for k and epsilon using Reynolds number.


Thank you

[mrishi]

1. If the basic flow characteristics are similar (similar velocity, pressure drop per unit length etc), then geometry should not change things very drastically, except in regions of curvature etc. I cannot really say at this point whether it is due to a poor mesh, or boundary conditions need to be set better. Have you confirmed that the new geometry is meshed well? (It is anyway good practice to run checkMesh on your case)

You can make decision to switch to turbulence based on value of Reynolds number for your pipe cross-section.


2. From the file description it says it should be able to handle turbulence. (https://www.openfoam.com/documentati...8C_source.html)

 Description
     Steady-state solver for buoyant, turbulent fluid flow and solid heat

     conduction with conjugate heat transfer between solid and fluid regions.

3. Look for solutions of flow in fully developed pipe flow. This will give you an idea of flow characteristics. You would find expressions such as:

k = 3/2(UI)^2
I = I(Re)
epsilon = C_mu( k^3/2)/l

These are sort of empirical correlations for turbulent flow in simple geometries. These values will give you realistic estimates of flow at boundaries, so that your solution in the domain stays realistic.


In your simplified heat exchanger, was your rectangular pipe also interacting with the heat sources? It seems to be located very differently for the purpose of heat exchange.

[Raza Javed]

Thank you for your prompt reply.


1. Yes, I have RUN the checkMesh utility to check my mesh. Here is the log for that:

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

Now, I am trying to calculate the Reynolds number of this new geometry.
Here I have one question that, while calculating the Reynolds number, what would be the characteristic length? Is it the full length of pipe(from inlet to outlet) OR is it just the diameter of the pipe?


Because if I take the full length of pipe then the Reynolds number turns out to be Turbulent, but If I select just the diameter of the pipe as a characteristic length, then it says the laminar.



In my simplified geometry the pipe was not touching the heat sources, but there is a small gap in between.

[mrishi]

There are quite a few elements with non-orthogonal faces, so that may be one of the areas to look into for potential improvement.


The diameter would be characteristic length so as to get an idea of flow through the cross-section. That is why I said changing geometry would not change things by much. Basically, I suppose you dont have to worry about turbulence at this point. The problem with your setup seems outside of this aspect.

Well, one thing you could try is as follows: in changing the case you have changed the geometry of pipe and its heat exchange characteristics with the plates.

Because you mention that your case setup works for the simple geometry (even though we may have scope to improve those BCs), it makes me think that if the problem comes because of temperature exchange at the fluid-solid contact and not because of the new mesh, then, if you used the simple rectangular pipe geometry, but brought it into contact with the walls, (similar to your actual case), then the simple geometry will also start to show similar error. This would mean that prime source of error is because the heat exchange interfaces are ill-defined. Basically this way we can isolate the effect of new mesh, non-orthogonality, change in length etc and instead first perfect the simpler model which is working on some basic level, then it would be easier to take the setup to the other geometry. I would suggest to investigate that pipe geometry with similar contact as the final geometry, then you could perhaps see whether the error appears on change of pipe geometry, or due to the kind of heat exchange setup. If the original setup still runs to steady state, then it would mean that it is the new geometry which is causing divergence. In this case we may want to start looking at the non-orthogonal cells and try nNonOrthogonalCorrectors 1 to improve results.

But move in small well-defined step changes that is very important. Do not make too many changes at once. It took me about a week or so to get my first simulation to become run, and that involved tweaking all these aspects.

[Raza Javed]

Originally Posted by mrishi

There are quite a few elements with non-orthogonal faces, so that may be one of the areas to look into for potential improvement.


The diameter would be characteristic length so as to get an idea of flow through the cross-section. That is why I said changing geometry would not change things by much. Basically, I suppose you dont have to worry about turbulence at this point. The problem with your setup seems outside of this aspect.

Well, one thing you could try is as follows: in changing the case you have changed the geometry of pipe and its heat exchange characteristics with the plates.

Because you mention that your case setup works for the simple geometry (even though we may have scope to improve those BCs), it makes me think that if the problem comes because of temperature exchange at the fluid-solid contact and not because of the new mesh, then, if you used the simple rectangular pipe geometry, but brought it into contact with the walls, (similar to your actual case), then the simple geometry will also start to show similar error. This would mean that prime source of error is because the heat exchange interfaces are ill-defined. Basically this way we can isolate the effect of new mesh, non-orthogonality, change in length etc and instead first perfect the simpler model which is working on some basic level, then it would be easier to take the setup to the other geometry. I would suggest to investigate that pipe geometry with similar contact as the final geometry, then you could perhaps see whether the error appears on change of pipe geometry, or due to the kind of heat exchange setup. If the original setup still runs to steady state, then it would mean that it is the new geometry which is causing divergence. In this case we may want to start looking at the non-orthogonal cells and try nNonOrthogonalCorrectors 1 to improve results.

But move in small well-defined step changes that is very important. Do not make too many changes at once. It took me about a week or so to get my first simulation to become run, and that involved tweaking all these aspects.

Hi...


Following your suggestion, I made a new simple geometry, and it is running without errors, and as I increase the simulation time, it is converging to steady-state.


I have attached picture of trials I have done at different simulation times (10s, 20s and 30s). you can see the temperature change in the fluid region due to heat sources, and After certain time, the whole fluid region and the plates take the same temperature.


For your reference, I am attaching my boundary conditions here also:

boundary
{
    inlet
    {
        type            patch;
    }
    outlet
    {
        type            patch;
    }
}
U
{
    internalField   uniform (0 0 1e-3);

    boundaryField
    {
        inlet
        {
            type                fixedValue;//pressureInletVelocity;//
            //volumetricFlowRate  0.066;
            //extrapolateProfile  yes;
            value               uniform (0 0 0.2);
        }

        outlet
        {
            type            zeroGradient;
            //value           uniform (0 0 0);//$internalField;
        }

       /* ".*"
        {
            type            noSlip;
        } */
        "fluid_to_box"
        {
            type            noSlip;
        }
    }
}

T
{
    internalField   uniform 300;

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

        outlet
        {
            type            zeroGradient;//inletOutlet;
            value           $internalField;
            //inletValue      $internalField;
        }

        /*".*"
        {
            type            zeroGradient;
            //value           uniform 300;
        }*/

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


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
    {
        inlet
        {
            type            zeroGradient;//fixedValue;
            //value           uniform 10;
        }
        outlet
        {
            type            fixedValue;
            value           uniform 0;
        }
        "fluid_to_.*"
        {
            type            zeroGradient;

        }
    }
}

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

[mrishi]

I am assuming this is a laminar solution as we have established that.
What I find strange is: is it normal to see the temperature of heat sources change as well?
Anyway, so in this exact situation keeping everything else constant, if instead of the central rectangular pipe, we switch to a) cylindrical pipe (straight), and b) cylindrical pipe (S shaped curved) they all should technically run fine, as long as your direction of velocity flow is correctly defined (it should be moving into the inlet face - might be important to consider when using curved pipe). What do you think?

[Raza Javed]

Originally Posted by mrishi

I am assuming this is a laminar solution as we have established that.
What I find strange is: is it normal to see the temperature of heat sources change as well?
Anyway, so in this exact situation keeping everything else constant, if instead of the central rectangular pipe, we switch to a) cylindrical pipe (straight), and b) cylindrical pipe (S shaped curved) they all should technically run fine, as long as your direction of velocity flow is correctly defined (it should be moving into the inlet face - might be important to consider when using curved pipe). What do you think?

Hi...


As you suggested, I made one geometry with cylinder (image attached). and it is working fine.


Then I changed the shape of cylindrical pipe. and then first it took a bit longer time in mesh generation. After that, running of simulation was also very slow, and then I got the following results.(image attached)


Now, what do you think about these results?


And one additional question here is that, how we check the direction of velocity?`By orientation of the geometry in Cartesian coordinate? for example, pipe is located in x direction, inlet is at (x=0) and outlet is at (x=10), then my velocity will be like this:


type fixedValue
value uniform (3 0 0);


is it like this OR something else?

[mrishi]

how we check the direction of velocity?`By orientation of the geometry in Cartesian coordinate?

Yes in global cartesian coordinates. If we assume the upstream flow is fully developed then it would enter normal to the inlet plane (the standard laminar pipe flow).

[Raza Javed]

Hi...


I think the problem occurs when we try to model turbulent flow, because in curved shape pipe, I think it doesn't remain laminar anymore, and we modeled it as a laminar flow, that is why we are unable to see the expected results.


OR do you think there could be some other problem?



Now, the problem is how can we model turbulent flow in this curved shape pipe, and what are basics? Here is my understanding, please correct me if I am wrong:


1. We need to define a certain inlet velocity so that we can find the Reynolds number?
2. Once we have the Reynolds number, then we need to calculate the boundary conditions for k and epsilon?
3. And after that, we need to follow the steps that you mentioned in earlier posts relating to turbulent flows?


kindly guide me, if I am missing something.


Thank you

[mrishi]

Why is the temperature in your curved pipe concentrated so asymmetrically? Have you double checked your direction of flow is correctly aligned in the new geometry? It is hard to tell based on your images. Can you check velocity vectors around inlet that it is oriented correctly?



If all that is correct then your curved pipe should also give a profile similar to the straight one if you have setup your inlet and outlet correctly and your new mesh is alright.



I feel that if your velocity and pipe diameter are not changing much then neither should the Reynolds number, but you can/should confirm this about curved pipes, I do not know much about that aspect but it should be available in pipe and aerodynamics studies etc.

[Raza Javed]

Originally Posted by mrishi

Why is the temperature in your curved pipe concentrated so asymmetrically? Have you double checked your direction of flow is correctly aligned in the new geometry? It is hard to tell based on your images. Can you check velocity vectors around inlet that it is oriented correctly?



If all that is correct then your curved pipe should also give a profile similar to the straight one if you have setup your inlet and outlet correctly and your new mesh is alright.



I feel that if your velocity and pipe diameter are not changing much then neither should the Reynolds number, but you can/should confirm this about curved pipes, I do not know much about that aspect but it should be available in pipe and aerodynamics studies etc.

Just an additional question here related to the flow direction, as you can see in the figure attached, the straight side of pipe is inlet and the curved side is outlet. And in the figure, you can see the Cartesian axis on the left bottom side. So, from this, can we say that our velocity at inlet is in the positive x-direction? that we can define like this:


value uniform (velocity 0 0);

[mrishi]

value uniform (velocity 0 0);

Yes this seems absolutely correct. If the case ran on the straight pipe it should ideally run on this curved pipe as well, because it is only an extension of the original shape, nothing else.



With this case, are you still getting the same error as your initial case, asked in first post? Otherwise it is quite strange if you actually got a solution with 1000 degrees at inlet and zero everywhere else, it seems very unphysical. In that case, the error had persisted with and without turbulence model, which is why I was concerned about the quality of BC and mesh setup.


If your error is still there, it seems that either there is some flaw introduced due to the new mesh with curvature related issues in quality (remember 4000+ non-orthogonal faces?), or the physics of the case have changed due to tendencies of curvature in terms of flow and temperature.


At this point, I am sorry that I can't precisely tell you what is happening, I am also curious to know. I am pretty sure it isn't a horrendous mistake we are making but a rather simple issue that needs to be fixed. See if playing with the mesh gets you somewhere. Also read about outerCorrectors and nonOrthogonalCorrectors they might be helpful in stabilizing your solution. Wherever in doubt, try replacing the BCs you don't fully understand with simple ones like fixedValue and zeroGradient. These initial divergences are typically mesh/BC problems.

[Raza Javed]

Thank you so much for your suggestion.


I will try to re-generate my mesh.



In the meanwhile, I have one question related to Reynold's number.


As we know the relation for Reynold's number:


Re=(rho*u*L)/v


where



rho = density of the fluid (1000 in my case)

u = velocity of fluid (1e-3 in my case)

L=Diameter of pipe(2.5mm in my case)
v= kinematic Viscosity (959e-6 in my case)



Using the above values, my Reynolds number turns out to be less than 2300, it means that my flow will be Laminar.


But if due to some reasons, the velocity of fluid changes somewhere in the pipe, would it change the fluid flow from Laminar to Turbulent? (because as per calculation using the above relation, if velocity increases more than 1.76m/s in my case, then the Re >4000, which means it is Turbulent flow)

[mrishi]

Technically you are right, it will be a low turbulence situation. However, velocity going from 1e-3 overall to 1.76 m/s in a local region isn't very common, that is a huge acceleration if you think about it - may occur if you change the cross-section suddenly by a large amount.



It may be possible if your average flow velocity itself is higher, then you'll have to use a turbulence model because your effective viscosity will have significant contribution of turbulence effect in addition to the normal laminar kinematic viscosity (959e-6). This will change the flow patterns substantially and may even be necessary for convergence.

[cryabroad]

Hi Raza,

I think your problem is more related to the physics, not modelling (simulation).

For example, is the fluid incompressible or compressible? If compressible, I assume that you are using some kind of equation of state to get density from temperature. I ask this question because in your output file, the maximum density reaches maxRho (2.0) really fast. Why did you bound your densities between 0.2 and 2.0? If the actual density of your fluid should be, say 5.0 (from some equation of state), but you are limiting it to be no larger than 2.0, of course your simulation blows up.

Second, I myself is not familiar with chtMultiRegionSimpleFoam, but based on the description it is a steady-state solver. But, in your output file you are having Time = 0.1, 0.2, 0.3, which means you are not using deltaT = 1 in your controlDict. In this case of course your solver blows up, because you are using it to solve a transient case, which it is not designed for.

Third, regarding the turbulence issues. The first thing you should do is to get a physical grasp of the type of the flow, mainly based on the Reynolds number. The curved geometry should not change the flow type, given your cross section area doesn't change. Curved geometry does give you a secondary flow, but that's not related to whether the flow being laminar or turbulent. With the geometry out of the way, here comes the real question-is your fluid incomressible or not? If it's compressible, its density changes with temperature, which means the density changes while heat transfer is happening. In this case, the density in the heated (or cooled) section might be very very different from what you have at the inlet. As a result, the flow changes from laminar to transition, and may eventually to turbulence. I don't really think this is what your case is about, but it's something to think about.

Again, you should provide more information regarding the physics of your problem, the modelling or simulation part is secondary. I suggest the following: check if the fluid is incompressible or not. if it's incompressible, it would be much easier and we go from there; if it's compressible, what equation of state best describe your fluid, what about other thermophysical properties, such as viscosity?

I feel like I asked a lot of questions, but hey, that's what you should do in CFD. Don't get lost in the modelling part of your problem (people do that a lot), start from fundamentals and physics. Figuring out how to condition your solver to match the physics is easy (although not as easy as you think), but figuring out the physics is much harder.

Regards,
Ruiyan

[Raza Javed]

Hi Ruiyan,


Thank you so much for your reply.


I have incompressible fluid in my simulation.


And I have changed my controlDict file, you can see below:

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

application     chtMultiRegionSimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         100;

deltaT          1;

writeControl    timeStep;

writeInterval   1;

purgeWrite      5;

writeFormat     ascii;

writePrecision  7;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

But after this, my simulation stopped at T=35s.


And the error is below:

Time = 35


Solving for fluid region fluid
DILUPBiCG:  Solving for Ux, Initial residual = 0.6843924, Final residual = 0.02674027, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.6412737, Final residual = 0.03455567, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.5605676, Final residual = 0.02308348, No Iterations 2
DILUPBiCG:  Solving for h, Initial residual = 0.6169205, Final residual = 0.01893477, No Iterations 3


--> 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)

I am also attaching the log file for your reference.