[bullmut]

Hi Hi


So i am trying to simulate a engine cell.
I am not looking at the engine itself but the exhaust gas.
So i have an inlet where the air gets drawn in, and outlet where the hot gas is expelled.
The room has a single inlet/outlet (called a port). The rest of the boundary conditions are walls.
I have attached an image for reference.


The problem is that teh denisty goes negative after a few iteratiorn causing a crash.
The following are the boundary condistions being used.

U boundaryField
{
    walls
    {
        type            noSlip;
    }
    inlet
    {
        type                flowRateOutletVelocity;
        volumetricFlowRate  0.49;
        value               uniform (0 0 0);
    }
    outlet
    {
        type                flowRateInletVelocity;
        volumetricFlowRate  1.11;
        value               uniform (0 0 0);
    }
    port
    {
    type zeroGradient;    
    }    
}

T - boundaryField
{
walls
    {
        type            zeroGradient;
    }
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 673;
    }
    port
    {
        type            inletOutlet;
        inletValue      uniform 300;
        value           uniform 300;
    }
}

p_rgh - boundaryField
{

walls

{
        type            fixedFluxPressure;
        gradient        uniform 0;
        value           uniform 101325;
    }
inlet & outlet
    {
        type            fixedFluxPressure;
        gradient        uniform 0;
        value           uniform 101325;
    }
port
    {
        type            fixedValue;
        value           uniform 101325;
    }
}

p is calculated from p_rgh


I know i have problem with my boundary conditions but i don't know hat they should be... any suggestions

[Bloerb]

Why are you using flowRateInletVelocity and flowRateOutletVelocity. And btw why outlet boundary conditions on the inlet and outlet boundary conditions at the outlet? Setting the velocity at inlet and outlet is pretty unstable. Typically you set the velocity at the inlet and use one of these at the outlet:

zeroGradient <- might cause problems with backflow
inletOutlet <- explicitly specifing the backflow (to e.g zero)
pressureInletOutletVelocity <- more suited if backflow can occur

The main reason for why this is unstable is the behaviour towards numerical error. If you specify velocity at inlet and outlet as boundary conditions your solver must find a solution which produces exactly these values. If you only set one of them, errors can be "flushed out" at the outlet and won't build up. Setting the velocity at the outlet is also possible but again you can't flush the errors out as easily since they are transported downstream.
You can specify boundary conditions like that, you can also use pressure pressure boundary conditions, but we need to know why you are using this combination.


Because you can try two ways to solve this:

1. improve your mesh and use schemes that are more stable
2. use boundary conditions that provide a higher stability

[Tobermory]

Bloerb - I don't think you understand his geometry. Have a look at the figure - the inlet and outlet are boundaries on the engine that is inside the room. The port boundary is the actual domain inlet/outlet.

Bullmut - where does the density go negative? I would suggest looking at where the density goes negative and then look at the p, U, T field values there to try work out which field is causing the trouble.

[bullmut]

@Bloerb, thanks for your reply, Ill try and address each of your concerns:

Why are you using flowRateInletVelocity and flowRateOutletVelocity.

I was given the flow rate specifications for the engine.

And btw why outlet boundary conditions on the inlet and outlet boundary conditions at the outlet

The engine outlet is a "inlet" to the simulation and the inlet is an "outlet" relative to the simulation. As Tobermory pointed out the only true outlet its called a port. I will try setting a pressure BC on the engine outlet (instead of a velocity one) and see if it improves stability.


@Tobermory, i actually tried to find this out. but i dont have the density file in my folders. I only see the error in the log.txt (see below)

SIMPLE: no convergence criteria found. Calculations will run for 1000 steps.

Reading thermophysical properties

Selecting thermodynamics package
{
type heRhoThermo;
mixture pureMixture;
transport const;
thermo hConst;
equationOfState perfectGas;
specie specie;
energy sensibleEnthalpy;
}

Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
RAS
{
RASModel kEpsilon;
turbulence on;
printCoeffs on;
Cmu 0.09;
C1 1.44;
C2 1.92;
C3 0;
sigmak 1;
sigmaEps 1.3;
}


Reading g

Reading hRef
Calculating field g.h

Reading field p_rgh

No MRF models present

Radiation model not active: radiationProperties not found
Selecting radiationModel none
No finite volume options present

Starting time loop

Time = 1

DILUPBiCGStab: Solving for Ux, Initial residual = 1, Final residual = 0.000206291, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 1, Final residual = 0.000248382, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 1, Final residual = 0.00901985, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 1, Final residual = 0.00225547, No Iterations 2
GAMG: Solving for p_rgh, Initial residual = 0.987749, Final residual = 0.00866302, No Iterations 5
GAMG: Solving for p_rgh, Initial residual = 0.00707345, Final residual = 4.85276e-05, No Iterations 9
GAMG: Solving for p_rgh, Initial residual = 4.83479e-05, Final residual = 3.23256e-07, No Iterations 10
time step continuity errors : sum local = 1.16471e-08, global = 1.37673e-09, cumulative = 1.37673e-09
rho max/min : 1.17641 0.524403
DILUPBiCGStab: Solving for epsilon, Initial residual = 0.0277441, Final residual = 6.02074e-06, No Iterations 1
DILUPBiCGStab: Solving for k, Initial residual = 1, Final residual = 0.000258576, No Iterations 1
ExecutionTime = 17.5 s ClockTime = 18 s

Time = 2

DILUPBiCGStab: Solving for Ux, Initial residual = 0.238591, Final residual = 4.62849e-05, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.261385, Final residual = 6.29764e-05, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.256434, Final residual = 6.25529e-05, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 0.105216, Final residual = 0.00018284, No Iterations 2
GAMG: Solving for p_rgh, Initial residual = 0.919036, Final residual = 0.00571887, No Iterations 5
GAMG: Solving for p_rgh, Initial residual = 0.0107584, Final residual = 5.73898e-05, No Iterations 12
GAMG: Solving for p_rgh, Initial residual = 5.75714e-05, Final residual = 4.72643e-07, No Iterations 10
time step continuity errors : sum local = 6.51067e-08, global = -1.93338e-09, cumulative = -5.56651e-10
rho max/min : 1.25484 0.00612429
DILUPBiCGStab: Solving for epsilon, Initial residual = 0.00161202, Final residual = 5.51097e-07, No Iterations 1
DILUPBiCGStab: Solving for k, Initial residual = 0.281582, Final residual = 8.23178e-05, No Iterations 1
ExecutionTime = 30.59 s ClockTime = 31 s

Time = 3

DILUPBiCGStab: Solving for Ux, Initial residual = 0.146195, Final residual = 2.56682e-05, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.113337, Final residual = 2.51388e-05, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.128397, Final residual = 2.82904e-05, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 0.0768671, Final residual = 0.000121801, No Iterations 2
GAMG: Solving for p_rgh, Initial residual = 0.727446, Final residual = 0.00701196, No Iterations 10
GAMG: Solving for p_rgh, Initial residual = 0.00444474, Final residual = 2.38272e-05, No Iterations 14
GAMG: Solving for p_rgh, Initial residual = 2.39488e-05, Final residual = 2.0509e-07, No Iterations 10
time step continuity errors : sum local = 4.24781e-08, global = 1.28846e-09, cumulative = 7.31805e-10
rho max/min : 1.20518 -0.574073
DILUPBiCGStab: Solving for epsilon, Initial residual = 0.00239489, Final residual = 7.70594e-07, No Iterations 1
bounding epsilon, min: -1.10268 max: 3013.88 average: 11.5738
DILUPBiCGStab: Solving for k, Initial residual = 0.1475, Final residual = 3.72107e-05, No Iterations 1
bounding k, min: -0.446489 max: 79.3595 average: 0.961464
ExecutionTime = 45.28 s ClockTime = 46 s

Time = 4

DILUPBiCGStab: Solving for Ux, Initial residual = 0.236749, Final residual = 6.89528e-16, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.276004, Final residual = 3.72831e-14, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.0854416, Final residual = 7.3865e-15, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 0.0594038, Final residual = 9.49646e-05, No Iterations 2

Since the density is negative the solver crashes during the p_rgh calculation.
How do i go about viewing the density in the paraview? Can i explicitly ask for rho to be saved in a time folder?


Thanks for the fast replies by the way!

[Tobermory]

You can force the density field (or other hidden fields) to be written using

buoyantSimpleFoam -postProcess -func 'writeObjects(rho)'

(Replace rho with bananas to see what fields are available - the solver will crash but will then list the acceptable options). Try getting the code to save after iter 3 and then investigate where it's falling over. Good luck!

[bullmut]

Originally Posted by Tobermory

You can force the density field (or other hidden fields) to be written using

buoyantSimpleFoam -postProcess -func 'writeObjects(rho)'

(Replace rho with bananas to see what fields are available - the solver will crash but will then list the acceptable options). Try getting the code to save after iter 3 and then investigate where it's falling over. Good luck!

Hehe, found this (How to get density-field for compressible flow?) and am doing just that! Will keep you posted
Cheers

[bullmut]

Ok i am back,



Negative density is reported at the engine inlet
Weirdly enough the timing is off. According to the log i get a negative density at step 3 but when i plot the values per time step i actually have a negative density at step 2.
I have not yet adjusted the BCs as suggested. Mainly because if i remove prescribed velocity BC and make it pressure.. what do i do? What should p_rgh be?
This is the situation as it stands

FYI - TURBINE_out == engine out & TURBINE_IN == engine inlet
when comparing code and images.

p_rgh

    "engine_outlet|engine_inlet"
    {
        type            fixedFluxPressure;
        gradient        uniform 0;
        value           uniform 101325;
    }

U

   engine_inlet
    {
        type                flowRateOutletVelocity;
        volumetricFlowRate  0.5;
        value               uniform (0 0 0);
    }
    engine_outlet
    {
        type                flowRateInletVelocity;
        volumetricFlowRate  0.5;
        value               uniform (0 0 0);
    }

[Paul Pickering]

I am not sure what the solution is to this problem, but I am having a similar problem with negative densities in buoyantSimpleFoam. I have a fairly complicated geometry and there is a crevice in the surface where negative densities occur and crashthe simulation on an early iteration. The equation of state being used is the perfect gas equation.


I am wondering if there is a way to put minimum limits on density (and perhaps other parameters) to allow the simulation to proceed to the next iterate? Does anyone know whether this is possible?

[Tobermory]

Originally Posted by bullmut

Negative density is reported at the engine inlet
Weirdly enough the timing is off. According to the log i get a negative density at step 3 but when i plot the values per time step i actually have a negative density at step 2.

Ok - that helps narrow down which BC is the culprit. From your log file, the min density on iter 1 matches your 673K engine outlet, but has become unphysical at iter 2. I wonder whether trying to fix the mass flow at that inlet rather than the volume flow might be easier for the code ... not as "correct" as a vol flow BC, though, since the engine will draw a constant vol flow. Have you tried playing with the underrelaxation factors - incl rho, to see if that helps? Alternatively, on some cases I have found that slowly ramping up a BC using the uniformFixedValue boundary can help the solver get started.

I believe the SIMPLE solver does struggle often if the initial field is far from the real solution, so if the above does not work, try initialise with potentialFoam & see if that helps.

I have not yet adjusted the BCs as suggested. Mainly because if i remove prescribed velocity BC and make it pressure.. what do i do? What should p_rgh be?

Indeed. You would have to dynamically change the boundary pressure value until you got the vol flow that you wanted ... possible I guess, but there HAS to be an easier way.

Best of luck - and let us know when you have cracked it!

[bullmut]

For those interested.. i have manged to get my simulation running.
the weird part is, it runs in openfoam 6 but not 8.. i think this is because the consolidated buoyantBoussinesqSimplefoam and buoyantSimpleFoam into 1 solver... in version 7 already.


so if you still want to know, i changed the boundary conditions as follows (paste from patchSummary):

Valid fields:
volScalarField nut
volVectorField U
volScalarField k
volScalarField p_rgh
volScalarField alphat
volScalarField p
volScalarField T
volScalarField epsilon

group : wall
scalar nut generic
scalar k generic
scalar p_rgh fixedFluxPressure
scalar alphat generic
scalar p calculated
scalar T zeroGradient
scalar epsilon generic
vector U noSlip

patch : TURBINE_IN
scalar nut calculated
scalar k zeroGradient
scalar p_rgh fixedFluxPressure
scalar alphat calculated
scalar p calculated
scalar T zeroGradient
scalar epsilon zeroGradient
vector U fixedValue

patch : TURBINE_OUT
scalar nut calculated
scalar k turbulentIntensityKineticEnergyInlet
scalar p_rgh fixedFluxPressure
scalar alphat calculated
scalar p calculated
scalar T fixedValue
scalar epsilon generic
vector U flowRateInletVelocity


patch : EXHAUST_OUT
scalar nut calculated
scalar k zeroGradient
scalar p_rgh fixedValue
scalar alphat calculated
scalar p calculated
scalar T zeroGradient
scalar epsilon zeroGradient
vector U zeroGradient

Time = 0

By setting the outlet (exhaust) to zeroGradient, prescribing the domain outlet velocity (turbine _in) and assigning a mass flow rate to the domain inlet (turbine_out) i have BC's that work.
The simulation runs for a few hundred iterations, then i get a spike nine the k and epsilon values. If i let it run on it eventually crashes. Not due a negative density though. rather a uncontrolled k or epsilon value. so the spike warns me runaway is coming.
I took the approach of turning turbulence off when i get the spike, and keeping the simulation going. 20 iterations later i turn on the turbulence again. I repeat this process 3 times before i leave turbulence on and run till i get convergence.



so i dont have the prefect answer - BUT i have a working solution that gives realistic results.



I am open to suggestions as to why my turbulence model spikes and causes a runaway situation.
Also any info on how the consolidated the solvers between 6 and 7.. i looked at github but i am not expert in reading code.


thanks for the suggestions everyone

[Tobermory]

Well done. I have found that OF can be really touchy about BCs and initial fields. For the problem with the ke spike, just check your initial fields incl nut. Do you have something credible (non-zero) set for nut? And credible k/eps values on the engine outlet? That has tripped be up before, for a simple jet simulation that kept on blowing up for no apparent reason ... almost drove me mad.