[Hr_kules]

Hey, recently i started using BuoyantSimpleFoam i did a small analysis for a pipe that worked exceptionally well, now i want to use another geometry. However if i run the case just with a different geometry (a corrugated channel instead of a pipe) the case is very instable and functionbojects as well as residuals are all over the place.
My boundary conditions are the following:

U

boundaryField
{

    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
    type fixedValue;
        value uniform (0.397441359254415 0 0);
    }

    outlet
    {
        type zeroGradient;
    }

    wallHeatedTop
    {
        type            noSlip;
    }
    wallHeatedBottom
    {
        type            noSlip;
    }    
    wallAdiabatic
    {
        type            noSlip;
    }   
}

T

boundaryField
{

    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
        type            fixedValue;
        value           uniform 298.15;
    }

    outlet
    {
        type            zeroGradient;
    }
    wallHeatedTop
    {
    type        fixedGradient;
    gradient    uniform 823.723229;
    }
    wallHeatedBottom
    {
    type        fixedGradient;
    gradient    uniform 823.723229;
    }    
    wallAdiabatic
    {
    type            zeroGradient;
    }    
}

p

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            calculated;
        value           $internalField;
    }
}

p_rgh

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type            fixedFluxPressure;
        value           $internalField;
    }
    outlet
    {
        type            fixedValue;
        value           $internalField;
    }    
}

omega

boundaryField
{

    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
    type         turbulentMixingLengthFrequencyInlet;
        mixingLength    0.000314802028337;
        value           $internalField;
    }
    outlet
    {
        type            inletOutlet;
        value           $internalField;
        inletValue      $internalField;
    }
    wallHeatedTop
    {
        type            omegaWallFunction;
        value           $internalField;
    }  
    wallHeatedBottom
    {
        type            omegaWallFunction;
        value           $internalField;
    }        
    wallAdiabatic
    {
        type            omegaWallFunction;
        value           $internalField;
    } 

}

k

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.061871583783885;
        value           $internalField;
    }
    outlet
    {
        type            inletOutlet;
        value           $internalField;
        inletValue      $internalField;
    }
    wallHeatedTop
    {
        type            kqRWallFunction;
        value           $internalField;
    }
    wallHeatedBottom
    {
        type            kqRWallFunction;
        value           $internalField;
    }          
    wallAdiabatic
    {
        type            kqRWallFunction;
        value           $internalField;
    }  
}

alpha

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
        type            calculated;
        value           $internalField;
    }

    outlet
    {
        type            calculated;
        value           $internalField;
    }
    wallHeatedTop
    {
        type            compressible::alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }
    wallHeatedBottom
    {
        type            compressible::alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }    
    wallAdiabatic
    {
        type            compressible::alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }     
    
}

nut

boundaryField
{

    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
        type            calculated;
        value           $internalField;
    }
    outlet
    {
        type            calculated;
        value           $internalField;
    }
    wallHeatedTop
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }
    wallHeatedBottom
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }    
    wallAdiabatic
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }
}

With following schemes and solution files:
fvSchemes

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)          bounded Gauss limitedLinear 1.0;
    div(phi,h)          bounded Gauss limitedLinear 1.0;
    div(phi,epsilon)     bounded Gauss limitedLinear 1.0;
    div(phi,omega)     bounded Gauss limitedLinear 1.0;
    div(phi,k)         bounded Gauss limitedLinear 1.0;
    div(phi,K)          bounded Gauss limitedLinear 1.0;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected 0.5;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected 0.5;
}
wallDist
{
    method meshWave;
}

fvSolution

solvers
{
    p_rgh
    {
        solver           GAMG;
        smoother         symGaussSeidel;
        tolerance        1e-7;
        relTol           0.01;
    }

    "(U|h|k|omega).*"
    {
        solver           PBiCGStab;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }
}

SIMPLE
{

    momentumPredictor yes;
    nNonOrthogonalCorrectors 6;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        U 1e-5;
        p 1e-4;
    }
    convergenceCriterion   1e-5;
    nNonOrthogonalCorrectors 0;  
}

relaxationFactors
{
    fields
    {
        rho             1.0;
        p_rgh           0.7;
    }
    equations
    {
        U               0.6;
        h               0.6;
        k               0.7;
        omega           0.7;
    }
}

Additionally a snippet of the misery:

DILUPBiCGStab:  Solving for Ux, Initial residual = 0.664815, Final residual = 0.0279068, No Iterations 1
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.609525, Final residual = 0.0235257, No Iterations 1
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.576135, Final residual = 0.0201185, No Iterations 1
DILUPBiCGStab:  Solving for h, Initial residual = 1, Final residual = 0.0342401, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 0.980315, Final residual = 2.79509e+45, No Iterations 1000
time step continuity errors : sum local = 2.61911e+76, global = 1.10847e+75, cumulative = 1.10847e+75
DILUPBiCGStab:  Solving for omega, Initial residual = 0.0335261, Final residual = 0.00216062, No Iterations 1
bounding omega, min: -8.37521e+26 max: 2.35929e+28 average: 1.8426e+23
DILUPBiCGStab:  Solving for k, Initial residual = 0.0322659, Final residual = 0.00302361, No Iterations 1
bounding k, min: -4.26595e+07 max: 5.82216e+08 average: 14976.7
ExecutionTime = 387.68 s  ClockTime = 388 s

wallHeatFlux wallHeatFlux1 write:
    writing object wallHeatFlux
    min, max, Q [W], q [W/m^2] for patch wallHeatedTop = 503.753, 503.753, 9.18111, 503.753
    min, max, Q [W], q [W/m^2] for patch wallHeatedBottom = 503.753, 503.753, 9.1811, 503.753

surfaceFieldValue t_fluid write:
    areaAverage(outlet) of T = 300.039

surfaceFieldValue outletTemperature write:
    areaAverage(outlet) of T = 299.953

surfaceFieldValue wallHeatedTopTemperature write:
    areaAverage(wallHeatedTop) of T = 297.311

surfaceFieldValue wallHeatedBottomTemperature write:
    areaAverage(wallHeatedBottom) of T = 297.341

surfaceFieldValue wallAdiabaticTemperature write:
    areaAverage(wallAdiabatic) of T = 297.232

yPlus yPlus1 write:
    writing object yPlus
    patch wallHeatedTop y+ : min = 2.85538e+35, max = 1.13867e+38, average = 3.53941e+36
    patch wallHeatedBottom y+ : min = 1.75043e+35, max = 1.46169e+38, average = 4.03137e+36
    patch wallAdiabatic y+ : min = 2.2173e+35, max = 8.39086e+37, average = 3.46597e+36

volFieldValue averageTemperatureFluid write:
    volAverage(channel) of T = 297.577

surfaceFieldValue deltaP.region1 write:
    areaAverage(inlet) of p = -1.64028e+93

surfaceFieldValue deltaP.region2 write:
    areaAverage(outlet) of p = 0

The checkMesh report stated no issues and the mesh was imported from gmsh.

I appreciate any hint or feedback!

[jenc24]

Hi,


did you try to run in laminar and see if the calculation runs smoothly?

[Cyrilub]

Hello Hr_kules,


I'm facing a similar issue, did you find any solution ?


Thanks

[piu58]

BuoyantSimpleFoam and its siblings are thought / tested for free convection. They accept addition external velocity, but have numerical problems with it. Normally, there is a small area which runs out of physics, and this spreads to the whole simulation.

I did not find any cure about this but another physical formulation.

[Hr_kules]

Hey, yes i did, i just changed the boundary for the velocity field to mass flow for some reason this worked in all the tested simulations and i got physically reasonable results.

Originally Posted by Cyrilub

Hello Hr_kules,


I'm facing a similar issue, did you find any solution ?


Thanks

[Hr_kules]

Originally Posted by piu58

BuoyantSimpleFoam and its siblings are thought / tested for free convection. They accept addition external velocity, but have numerical problems with it. Normally, there is a small area which runs out of physics, and this spreads to the whole simulation.

I did not find any cure about this but another physical formulation.

Yeah, since future plans involve more extensive simulations it is planned to work from simpleFoam and add desired pyhsical formulations.

[Cyrilub]

Already tried, it didn't worked for me, but thanks for the reply !