[FranckL]

Hello all,
It's my first post here, even if I've read the forum many times.
I'm specialized in pollution modeling, and it's been some weeks I have started to use CFD and especially openFoam.

I work on a problem and I'd like to have inputs from expert on good practices for my case. I want to model a gas leak from a pipe (methane) on an industrial site. I've tried many things, finally trying to derive the set-up from Fiates & Vianna 2016 ("Numerical Modelling of gas dispersion using openFoam").

I summarize what I did so far :
- I run rhoReactingBuoyantFoam
- fvScheme :

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    //default         none;

    div(phi,U)       Gauss upwind;
    div(phi,Yi_h)    Gauss upwind;
    div(phi,K)       Gauss upwind;
    div(phi,p)       Gauss upwind;
    div(phi,p_rgh)   Gauss upwind;
    div(phi,k)       Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega)   Gauss upwind;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected; 
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default no;
    p_rgh;
}

wallDist
{
    method          meshWave;
}

- fvSolution:

solvers
{
    "(rho|rhoFinal)"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
    };

    p_rgh
    {
        solver                  GAMG;
        tolerance               1e-5;
        relTol                  0.01;
        smoother                GaussSeidel;
        cacheAgglomeration      true;
        nCellsInCoarsestLevel   10;
        agglomerator            faceAreaPair;
        mergeLevels             1;
    };

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-6;
        relTol          0;
    };

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

    "(U|Yi|k|h|omega|epsilon)Final"
    {
        $U;
        tolerance         1e-6;
        relTol            0;
    };

}

PIMPLE
{
    momemtumPredictor           yes;
    nOuterCorrectors            1;
    nCorrectors                 2;
    nNonOrthogonalCorrectors    0;
}

- During the mesh generation, after blockMesh and snappyHexMesh, I add a "jet" bounday field with topoSet and createPatch that I will configure as an inlet (FYI, I ran checkMesh and it says it's OK):

Final polyMesh/boundary:
5
(
    walls
    {
        type            wall;
        inGroups        1(wall);
        nFaces          24582;
        startFace       2516067;
    }
    inletair
    {
        type            patch;
        nFaces          945;
        startFace       2540649;
    }
    outlet
    {
        type            inletOutlet;
        nFaces          6730;
        startFace       2541594;
        faces           ( ( 1 2 6 5 ) ( 4 5 6 7 ) ( 0 1 5 4 ) ( 3 2 6 7 ) );
    }
    building
    {
        type            wall;
        inGroups        1(wall);
        nFaces          91942;
        startFace       2548324;
    }
    jet
    {
        type            inlet;
        nFaces          4;
        startFace       2640266;
    }

- Zoom on the leak (in green):


- I set-up the boundary conditions for CH4, N2, O2, U...so that to simulate a leak after 25 s of run; example of CH4 file:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      CH4;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [ 0 0 0 0 0 0 0 ];

internalField   uniform 0;

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

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

    walls
    {
        type            zeroGradient;
    }

    building
    {
        type            zeroGradient;
    }

    jet
    {
        type             uniformFixedValue;
        uniformValue     table
        (
            (0.0   0)
            (25    0)
            (26    1)
            (625   1)
            (626   0)
        );
    }
}

Similar corresponding conditions for N2 and O2 so that N2+O2+CH4 always = 1.

Now my problems:
I tried different things for U as boundary conditions at the jet inlet:
a) Directly specifiying the flow speed
- If I use a high enough velocity (100 m/s) after 26s of run, i see the creation of the cloud and its further propagation on the site;
- If I limit the jet velocity to say 15/25 m/s (given the size of my leak on the mesh, it would correspond to a flow rate of a few hundreds L/min), the cloud remains very small, with very low CH4 concentration, a few cells around the leak only after 90s (I stopped the run here).
b) Directly specifying the flow rate with flowRateInletVelocity
- With the value I want to simulate (250 L/min, i.e. 4.2e-3 m3/s), same remark as above, almost no CH4 concentration is propagating.


Since I am an openFoam newbie, it's highly possible I have made basic errors...I'd be happy to hear from you if you see something wrong in my set-up, or if you have reference works/tutorials about gas leak dispersion with openfoam I have missed.


PS: I didn't show eveything (my post seems long to me already...), but I if you need more info on my set-up I'd be happy to share.

[Tobermory]

Can I clarify - what do you mean when you say

almost no CH4 concentration is propagating

?

Do you mean that CH4 concentrations are zero outside of the cells immediately surrounding the source (i.e. the CH4 is definitely not propagating)? If that is the case, check your solver output and the number of iterations that are being done on the CH4 equation. If they have dropped to zero, then adjust your solver tolerances to get the solver working again, or put a minIter line in the fvSolution file for CH4 etc.

However, if you mean there is a plume propagating, but that the concentrations are very low, then that's a whole different ballpark. The questions you need to ask youself are: is this correct behaviour (ie should it dilute so quickly)? Have I got enough mesh resolution to resolve the plume (coarse mesh resolution will smear your plume rapidly)? Are my numerical schemes overly diffusive (i.e. generating numerical or false diffusion ... you are using upwind, so the answer is yes)?

Hopefully that steers you in the right direction. Getting an appropriate mesh resolution of the pume can be quite a challenge in a "real 3D" scenario. Good luck!

[FranckL]

Thank you for your answer!

I indeed mean I only have small CH4 concentrations around the leak, but not zero. After a few second, I can reach 1 ppm in the surrounding cells, decreasing very rapidly.
My mesh resolution is increased in the jet area, but maybe not enough, so I will definitely check this. I may also have numerical diffusion...by the way, is it possible to check the intensity of numerical diffusion easily?

As to if it is a correct behavior...honestly I don't know yet. We have planned a field campaign, but not in the near future.

Also, I currently use the same setup for the whole run, but I thought I could first use a steady-state solver to reach an equilibrium in the domain, and then switch to the current set-up with rhoReactingBuoyantFoam and start the leak. This would save me time with my test, but I am not sure I can do it: is it possible to change solver during the run, with the mesh being the same obviously? I think if I start with the latest timestep it should be feasible. What do you think?

[Tobermory]

Ah - okay, then I suspect that focussing on the divScheme settings and mesh resolution is probably your best bet.

Upwind differencing (1st order) is nice and stable, but is very diffusive. Try some of the other higher order differencing schemes, but be aware - that may make your run slightly less stable. There's no easy way to identify the numerical error, since to do that you'd need the exact solution ... the simplest way though is to change a parameter and see if it makes a significant difference.

On the mesh - this is always the most challenging aspect. Again, the best way to prove to yourself that the mesh is not unduly affecting the solution is to refine the mesh and run again ... and keep going until the changes are negligible (or you run out of comuting resources!). Focus the refinement on the areas that you think are important to resolve accurately (e.g. the jet nozzle, the plume downstream), and add more refinement in there - you should start to see some real improvement. You can use topoSet to grab cells based on their concentration, and refineMesh to refine those cells and then mapFields to map the old solution onto the new mesh (so that you don't have to start the simulation from scratch).

Regarding the "correct solution", one of the things that I told my new CFD engineers is that it is always much easier to make a good simulation if you now what the answer should be! That's easy to say of course ... so how to estimate the right answer for your simulation? A field campaign is a great way, but won't help you now. Try looking for simple analogues in the literature (eg jet in a crossflow). You may only find data on velocity, but your plume velocity is a good measure of how your concentration is diluting (they both dilute at roughly the same rate), and so that is okay. Alternatively - look for a simpler plume model (EPA has some, for example; there are commercial models as well like PHAST .. expensive though). A quick run of one of these models will give you a feel for how quickly the plume should dilute.

Finally, can you restart a solution with a different solver? Sometimes yes ... providing that the new solver uses all the same variables as the old solver. If it requires new variables, eg p_rgh instead of p, then you may need to spend some time creating the starting files for these variables for the new solver.

[FranckL]

Thanks for your inputs, I will consider all of your advices to improve my simulation.