[ingraban]

I try to do a simulation of a fire in a tunnel with longitudinal air flow.
- OpenFoam Version v1812
- fireFoam case based on tutorial smallPoolFire3D

- The BCs at the inlet portal are
u: flowRateInletVelocity
p: calculated
p_rgh: prghTotalHydrostaticPressure
ph_rgh: fixedFluxPressure

- The BCs at the outlet portal are
u: pressureInletOutletVelocity
p: calculated
p_rhg: prghTotalHydrostaticPressure
ph_rgh: fixedValue

- The BCs at the tunnel walls are
u: noSlip
p: zeroGradient
p_rgh: fixedFluxPressure;
ph_rgh: fixedFluxPressure;

- The fire is defined similarly to the inlet portal, just giving a mass flow rate of CH4.

The setup works fine without the fire (giving a mass flow rate of air instead of CH4) and with a CH4 fire, but with zero flow rate at the inlet portal. I tried multiple variations of simple and derived boundary conditions for pressure and velocity.

With the fire and the air flow, the fire appears to affect the inlet condition. The longitudinal flow is not as defined in the BC. The visualisation of the flow velocity shows peaks close to the BC. The picture shows a longitudinal section though the domain, close to the inlet BC. We are looking at Ux (should be evenly distributed across the inlet). What am I doing wrong?

[ingraban]

I am still struggling with these boundary conditions. So far, I limited the fire size by reducing the area of the burner (but it is not as specified by the inlet flow of methane). And the longitudinal flow Ux is way above the specified flow.
Is it lack of convergence or an error in the boundary conditions? And what else could I try? Maybe the grid is too coarse? So far, I am limited with my computer power (~1 Mio. cells) and I don't want to tap external resources until I am sure about the simulation.

Courant Number mean: 0.2937 max: 0.605189
deltaT = 0.00833333
Time = 299.992

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
DILUPBiCGStab:  Solving for Ux, Initial residual = 7.70691e-06, Final residual = 3.64904e-08, No Iterations 1
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.000121998, Final residual = 5.9943e-07, No Iterations 1
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.000146818, Final residual = 8.45358e-07, No Iterations 1
DILUPBiCGStab:  Solving for O2, Initial residual = 8.69083e-05, Final residual = 6.61415e-07, No Iterations 1
DILUPBiCGStab:  Solving for H2O, Initial residual = 8.79887e-05, Final residual = 6.69055e-07, No Iterations 1
DILUPBiCGStab:  Solving for CH4, Initial residual = 6.22256e-05, Final residual = 1.4532e-07, No Iterations 1
DILUPBiCGStab:  Solving for CO2, Initial residual = 8.79887e-05, Final residual = 6.69055e-07, No Iterations 1
DILUPBiCGStab:  Solving for h, Initial residual = 9.88587e-05, Final residual = 7.16077e-07, No Iterations 1
min/max(T) = 282.591, 1594.31
GAMG:  Solving for p_rgh, Initial residual = 0.001461, Final residual = 2.49336e-05, No Iterations 2
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 8.4612e-08, global = -9.65967e-10, cumulative = 0.000115376
GAMG:  Solving for p_rgh, Initial residual = 0.000213141, Final residual = 6.35611e-07, No Iterations 3
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.15665e-09, global = 1.16067e-11, cumulative = 0.000115376
DILUPBiCGStab:  Solving for k, Initial residual = 1.32612e-05, Final residual = 5.22765e-08, No Iterations 1
bounding k, min: 0 max: 53.9932 average: 2.84801
ExecutionTime = 41790.8 s  ClockTime = 42541 s

Courant Number mean: 0.293699 max: 0.605172
deltaT = 0.00833333
Time = 300

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
DILUPBiCGStab:  Solving for Ux, Initial residual = 7.52108e-06, Final residual = 3.56617e-08, No Iterations 1
DILUPBiCGStab:  Solving for Uy, Initial residual = 0.000117884, Final residual = 5.67554e-07, No Iterations 1
DILUPBiCGStab:  Solving for Uz, Initial residual = 0.000139513, Final residual = 7.86985e-07, No Iterations 1
DILUPBiCGStab:  Solving for O2, Initial residual = 8.9583e-05, Final residual = 6.52377e-07, No Iterations 1
DILUPBiCGStab:  Solving for H2O, Initial residual = 9.06789e-05, Final residual = 6.59924e-07, No Iterations 1
DILUPBiCGStab:  Solving for CH4, Initial residual = 6.69555e-05, Final residual = 1.44643e-07, No Iterations 1
DILUPBiCGStab:  Solving for CO2, Initial residual = 9.06789e-05, Final residual = 6.59924e-07, No Iterations 1
DILUPBiCGStab:  Solving for h, Initial residual = 0.000101999, Final residual = 7.08459e-07, No Iterations 1
min/max(T) = 282.591, 1594.46
GAMG:  Solving for p_rgh, Initial residual = 0.00140299, Final residual = 2.46454e-05, No Iterations 2
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 8.36311e-08, global = -2.9129e-10, cumulative = 0.000115376
GAMG:  Solving for p_rgh, Initial residual = 0.000213343, Final residual = 6.3965e-07, No Iterations 3
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.17012e-09, global = 1.75397e-11, cumulative = 0.000115376
DILUPBiCGStab:  Solving for k, Initial residual = 1.3248e-05, Final residual = 5.21712e-08, No Iterations 1
bounding k, min: 0 max: 53.9931 average: 2.84801
ExecutionTime = 41792 s  ClockTime = 42542 s

volFieldValue HRR write:
    volIntegrate(region0) of Qdot = 3.97495e+07

End

Can anyone give me a pointer, please?

[ingraban]

It appears that the boundary condition for hydrostatic pressure is overly restrictive and causes the inlet flow rate to change from the defined velocity. The following works better, but still, the heat release rate is not what I expected.

- The BCs at the inlet portal are
u: fixedValue
p: calculated
p_rgh: fixedFluxPressure
ph_rgh: fixedFluxPressure

- The BCs at the outlet portal are
u: pressureInletOutletVelocity
p: calculated
p_rhg: prghTotalHydrostaticPressure
ph_rgh: fixedValue

- The BCs at the tunnel walls are
u: noSlip
p: calculated
p_rgh: fixedFluxPressure
ph_rgh: fixedFluxPressure

- The fire is defined similarly to the inlet portal, just giving a mass flow rate of CH4.

[clapointe]

I've seen your other thread too -- it sounds like you have figured out a working set of boundary conditions but the heat release is a problem. If you are still using infinitely fast chemistry, you might try switching to an eddy dissipation combustion model. This will couple your chemistry to your turbulence model. There are various flavors available in the FM Global repository for fireFoam (google fireFoam dev) and porting it to the version of foam that you are using shouldn't be too difficult. Otherwise, another option is to try and tweak the infinitely fast model to produce more realistic results (by changing the C parameter, for example). This is not to say that there could be other factors -- e.g. turbulence modeling, numerical schemes, mesh -- that are also impacting your results.

Caelan

[ingraban]

Thanks a lot for your suggestion. Yes, I also have some reservations against an Infinitely Fast Chemistry fire. Sounds like each time step is a small detonation. I'll be looking into the Eddy Dissipation Model right away. Hopefully, I'll find a suitable tutorial case...

[clapointe]

The FM global git repo has a couple of cases you should be able to use.

Caelan

[ingraban]

A search on "FM global git repo" produced the following link, which - again - is very helpful indeed. Thank you!

https://github.com/fireFoam-dev

[yambanshee]

Originally Posted by ingraban

I try to do a simulation of a fire in a tunnel with longitudinal air flow.
- OpenFoam Version v1812
- fireFoam case based on tutorial smallPoolFire3D

- The BCs at the inlet portal are
u: flowRateInletVelocity
p: calculated
p_rgh: prghTotalHydrostaticPressure
ph_rgh: fixedFluxPressure

- The BCs at the outlet portal are
u: pressureInletOutletVelocity
p: calculated
p_rhg: prghTotalHydrostaticPressure
ph_rgh: fixedValue

- The BCs at the tunnel walls are
u: noSlip
p: zeroGradient
p_rgh: fixedFluxPressure;
ph_rgh: fixedFluxPressure;

- The fire is defined similarly to the inlet portal, just giving a mass flow rate of CH4.

The setup works fine without the fire (giving a mass flow rate of air instead of CH4) and with a CH4 fire, but with zero flow rate at the inlet portal. I tried multiple variations of simple and derived boundary conditions for pressure and velocity.

With the fire and the air flow, the fire appears to affect the inlet condition. The longitudinal flow is not as defined in the BC. The visualisation of the flow velocity shows peaks close to the BC. The picture shows a longitudinal section though the domain, close to the inlet BC. We are looking at Ux (should be evenly distributed across the inlet). What am I doing wrong?

For future note: In general it is bad practice to define both pressure and velocity at an inlent and outlet for a single inlet/outlet case. A more stable solution will always be had by setting the velocity and pressure as a fixed value at one inlet/outlet, and specifying it's gradient at the other. They do not need to both be fixed at one inlet/outlet.

For example for a pipe flow case:

at inlet: fixed U and zero gradient P
at outlet: fixed P and zero gradient U
or
at inlet: fixed U and fixed p
at outlet: zero gradient p and zero gradient U

By using more than 1 fixed value for either pressure or velocity, you are over constraining your solution. For example:

inlet:
U = 10m/s
p = 0 Pa

outlet:
U = zeroGradient
p = 20 Pa

We know that normal pipe flow will cause a pressure drop, but by specifying the pressure at both you have over constrained the solution and no valid solution exists, which will cause the simulation to diverge.

However, I'm glad you managed to sort this out by yourself!

[ingraban]

Thank you for your response, yambanshee. I know about over-constraining of a pipe flow system. What's new for me here is the three pressure boundary conditions p, p_rgh and ph_rgh. As I understand, it is

p = ph_rgh + rho*gh + pRef and p_rgh = p - rho*gh

So, as a boundary condition, you have to define two of them and leave the third as "calculated", right?

[clapointe]

In fireFoam's pEqn, we see that p and p_rgh (not ph_rgh) are related as :

p = p_rgh + rho*gh + pRef;

So p is computed from p_rgh, rho, gh, and pRef. Indeed, we never actually solve for p -- it just is calculated from other known (and partially solved for) quantities to keep the thermo consistent. This -- as far as I know -- is why we have calculated bc's for p.

Now, we don't see ph_rgh in that equation -- that is because it is not used to compute p or p_rgh -- or at least not directly. It is used if "hydrostaticInitialization" is chosen by the user. Then a hydrostatic pressure field -- ph_rgh -- is computed before the simulation really starts. This field is then used by the prghTotalHydrostaticPressure boundary condition for the p_rgh field. This bc was added somewhere around foam version 4 to better model entrainment for large domains if I remember correctly. Without even looking at the math, we can guess that it might not be too beneficial for confined cases like the one we are discussing (with not much variation in a background hydrostatic pressure).

Now, having scanned the bc's again -- p should always be "calculated" because it is computed from other quantities. But otherwise the other bc's should work -- setting bc's for (compressible) heat transfer problems is always difficult. At the outlet, U is still pretty much zero gradient -- pressureInletOutletVelocity "guesses" a velocity for any faces with backflow (as determined by pressure) but defaults to zero gradient. At the inflow, velocity is specified and -- from my understanding of fixedFlux -- the pressure gradient is modified to be consistent with the velocity. So it does not appear to be over defined. The heat transfer solver tutorials use similar setups I think.

Caelan