[Jelle Lagerweij]

Dear (hopefully knowledgable) openFoam friends,
I am close to running a case which I have been working on for the last 2 weeks. However, there is this last error which I can't seem to resolve. Instead of a normal error message, indicating nicely where my error is, openFoam just seems to crash when I try to run my case. All the preparations seem to be ok, at least, snappyHexmesh returns the shape I expect, with all the correct in and outlet patches. And when I run the same problem in simpleFoam, it seems to work.


I switched to reactingFoam as I need to add mixing of Hydrogen in air, and off course, lots of extra inputs are needed for that. But the errors I get are the following:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
Build  : 9-c96968adf6d8
Exec   : reactingFoam -parallel
Date   : Jul 21 2022
Time   : 13:07:05
Host   : "calculus"
PID    : 146356
I/O    : uncollated
Case   : /home/xintc/OpenFOAM/xintc-9/run/Mixing
nProcs : 6
Slaves : 
5
(
"calculus.146357"
"calculus.146358"
"calculus.146359"
"calculus.146360"
"calculus.146361"
)

Pstream initialised with:
    floatTransfer      : 0
    nProcsSimpleSum    : 0
    commsType          : nonBlocking
    polling iterations : 0
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


PIMPLE: No convergence criteria found


PIMPLE: Operating solver in transient mode with 1 outer corrector
PIMPLE: Operating solver in PISO mode


Using LTS
Reading thermophysical properties

Selecting thermodynamics package 
{
    type            hePsiThermo;
    mixture         multiComponentMixture;
    transport       sutherland;
    thermo          janaf;
    energy          sensibleEnthalpy;
    equationOfState perfectGas;
    specie          specie;
}

Reading field U

Reading/calculating face flux field phi

Creating turbulence model.

Selecting turbulence model type RAS
Selecting RAS turbulence model kEpsilon
bounding k, min: 0 max: 0.12 average: 0.12
kEpsilonCoeffs
{
    Prt             0.85;
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              0;
    sigmak          1;
    sigmaEps        1.3;
}

Creating thermophysical transport model

Selecting thermophysical transport type RAS
Selecting default RAS thermophysical transport model unityLewisEddyDiffusivity
Creating reaction model

Selecting combustion model none
Creating field dpdt

Creating field kinetic energy K

No MRF models present

No fvModels present
No fvConstraints present
[2] #0  [3] #0  [4] #0  Foam::error::printStack(Foam::Ostream&)[5] #0  Foam::error::printStack(Foam::Ostream&)[0] [1] #0  Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)Foam::error::printStack(Foam::Ostream&)#0  Foam::error::printStack(Foam::Ostream&) at ??:?
 at ??:?
 at ??:?
 at ??:?
 at ??:?
[1] #1  Foam::sigFpe::sigHandler(int)[0] #1  Foam::sigFpe::sigHandler(int)[2] #1  Foam::sigFpe::sigHandler(int)[3] #1  Foam::sigFpe::sigHandler(int) at ??:?
[4] #1  Foam::sigFpe::sigHandler(int)[5] #1  Foam::sigFpe::sigHandler(int) at ??:?
[1] #2  ? at ??:?
[3] #2  ? at ??:?
[2] #2  ? at ??:?
[4] #2  ? at ??:?
[0] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"

I can see that the error occurs when I initialise the time loop, as a similar (although 2D) tutorial file runs with the following log.

No fvModels present
No fvConstraints present


Starting time loop

Time scales min/max:
    Flow        = 0.0001, 0.0001
    Temperature = 4.494232837e+307, 4.494232837e+307
    Composition = 4.494232837e+307, 4.494232837e+307
    Overall     = 0.0001, 0.0001
Time = 10

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab:  Solving for Ux, Initial residual = 1, Final residual = 2.459369251e-05, No Iterations 1

Can someone explain to me what such an Ostream error means and why it on general occurs? And how is that related to sigHandler?


PS: Maybe it is an important aspect, I do use adjustTimeStep with a maxCo of 0.1. I had another case that ran for a certain amount of time steps before it crashed with a similar error.

[Jelle Lagerweij]

I decided to open start fiddling around in the solver. So I copied the reactingFoam solver to my own run location and renamed it (using the

https://openfoamwiki.net/index.php/How_to_add_temperature_to_icoFoam

guide).


I started to include print statements to see in which line my code broke, it seems to be the following statement where it crashes:

turbulence->validate();

What would be the exact cause of this? Again, note that it crashes in a way that indicates illegal handling of floats of integers (like division by 0 or roots of negative values).


What does this line actually do? With kind regards,
Jelle Lagerweij

[dlahaye]

Hard to say from the data you share.

I suspect some error in the input (mesh, system folder, constant folder, 0 folder or elsewhere).

What is the output of checkMesh?

Can you share info on your case?

Does e.g. pimpleFoam (no chemistry) run fine?

[Jelle Lagerweij]

I added my problem as a zip. (if you execute ./run in the terminal it runs the case as I set it up).


I also did a checkMesh and I ran my problem using pimpleFoam. In pimpleFoam it runs fine, but I do have one skewed face. I would expect that there is an error in the initial fields and/or boundaries of N2 or H2.


Note that the mesh is not really appropriate for this case (I believe that my cells at the wall are too many wall distances away from it). But for now, I am fine with that. My goal was to just have it working correctly once before I would start implementing the real detailed case. In a similar way, I think that the turbulent intensities at the inlets are also not that correct in their value.


I hope that you can help me further.
With kind regard
Jelle Lagerweij


PS: I forgot to adjust the control dictionary. I compiled my own reacting foam solver (with extra print statements), so the application line in the controlDict should be changed from my_reactingFoam to reactingFoam.

[dlahaye]

Apologies, I am unwilling to run your case.

I can provide advice.

Does reactingFoam run with chemistry turned off?

[Jelle Lagerweij]

Oh, I need to explain the case in a little bit more detail then. I need to check if a room in which hydrogen is produced is vented enough. At an inlet I have a fixed velocity (Fan) and at an outlet, I have a fixed pressure. At the fan, we only have air (I called it N2 to make my live easier).



Additionally, I have an assumed leak of pure hydrogen (another inlet), with some fixed velocity.


I want to check how fast the high concentration of hydrogen is mixed with the air and, specifically, if at the outlet of the room the concentration is low enough to be safely released (sub 3 or 4%).


Therefore, I have the reactions turned off:

combustionModel  none;

active    false;

noCombustionCoeffs
{
}

By changing my own solver, I found that the error occurs just before the real time loop is started my code crashes. (in the line

turbulence->validate();

). This is when the first computation is done on the created fields. So I expect that I made an error in my 0 folder files. As the model runs with the same 0 folder in pimpleFoam, I expect that I need to turn my attention to the initiations of the concentrations:


For N2 it is:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       volScalarField;
    location    "0";
    object      N2;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 1;

boundaryField
{
    Inlet1
    {
        type            fixedValue;
        value           uniform 1;
    }
    Inlet2
    {
        type            fixedValue;
        value           uniform 0;
    }
    
    Outlet
    {
        type            inletOutlet;
        inletValue      uniform 1;
        value           uniform 1;
    }
    
    //Outlet
    //{
    //    type            zeroGradient;
    //}
    
    Wall
    {
        type            zeroGradient;
    }
}


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

and for H2 it is:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       volScalarField;
    location    "0";
    object      H2;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    Inlet1
    {
        type            fixedValue;
        value           $internalField;
    }
    
    Inlet2
    {
        type            fixedValue;
        value           uniform 1;
    }
    
    Outlet
    {
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }
    
    //Outlet
    //{
    //    type            zeroGradient;
    //}
    
    Wall
    {
        type            zeroGradient;
    }
}


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

Maybe I made an obvious error somewhere in here.


With kind regards Jelle


PS: I definitely understand being unwilling to run my code. It is, first of all, unsafe from infections and besides that always a hassle to get right (specific openFoam versions, different computer hardware setup ect.).

[dlahaye]

Dear Jelle,

Thank you for your understand. Apologies for the delay in this reply.

I humbly suggest to try to run rhoPimpleFoam in a next step.

In using reactingFoam, we have repeatedly overlooked setting the boundary condition for the density at the inlet patch appropriately.

Once rhoPimpleFoam runs appropriately, we can take a next step.

Does this help? Kind wishes. Domenico.

[Jelle Lagerweij]

Dear Domenico,
It helps a lot. I got the task to make this simulation (turbulent two species flow simulation in a 40ft container) during the internship of my master. I only got a 6 core CPU and 64GB of RAM to make it work. I did do some work in Ansys Fluent before and wrote my own laminar axisymmetric CFD code in python, but I have never used openFoam before this. You seem to understand some parts of the inner working, building up the problem to an increasingly more complex case seems to be a much better approach than to start with reactingFoam immediately.

At least, we now found the issue at hand. I get the same type of error reporting, in the same place in the solver. While tweaking around for the correct thermophysicalProperties file for the rhoPimpleFoam case, I directly found how I could reduce some computational power in the future. I used Sutherland and Janaf coefficients before, however my case has no reactions, so no large temperature changes are expected. In the rhoPimpleFoam tutorial files, I found examples with const and hconst models for transport and thermodynamical properties respectively. Therefore, I rewrote my thermophysicalProperties file to the following:

FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            hePsiThermo;
//    mixture         multiComponentMixture;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    specie          specie;
    equationOfState perfectGas;
    energy          sensibleEnthalpy;
}

//species         2 (AIR H2);
//defaultSpecie AIR;

mixture
//AIR
{
    specie
    {
        molWeight   28.9;
    }
    thermodynamics
    {
        Cp          1005;
        Hf          0;
    }
    transport
    {
        mu          1.82e-05;
        Pr          0.71;
    }
}

//H2 // For now, H2 has same properties as Air of course this should be improved lateron
//{
//    specie
//    {
//        molWeight   28.9;
//    }
//    thermodynamics
//    {
//        Cp          1005;
//        Hf          0;
//    }
//    transport
//    {
//        mu          1.82e-05;
//        Pr          0.71;
//    }
//}


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

(note that in reactingFoam, I will use the multiComponentMixture, but I turned that, and the H2 properties off for now.)

For sofar, I understood that going from pimpleFoam to rhoPimpleFoam, I need to add the following additional boundary and initial conditions, for T and for alphat. As rhoPimpleFoam does not use multispicies transport yet, I will not have to describe the species distribution. I however get the same error, so it seems that my problem lies with these conditions (or the thermophysical properties file).

Lets first show the temperature file. I believe that this one should be correct. I would expect openFoam to compute the density as function of T and p as I used perfectGas as EOS.

FoamFile
{
    format      ascii;
    class       volScalarField;
    location    "0";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 292;

boundaryField
{
    Inlet1
    {
        type            fixedValue;
        value           $internalField;
    }
    
    Inlet2
    {
        type            fixedValue;
        value           $internalField;
    }
    
    Outlet
    {
        type            zeroGradient;
    }
    
    Wall
    {
        type            zeroGradient;
    }
}


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

The alphat initial and boundary conditions is the one I doubt most. I am not fully sure which characteristic is set by this function and what it does. I set these conditions in the same way as I found in some tutorials, without understanding (probably not the right approach). I believe it is connected to the heat transfer at (turbulent) wall.

I must note that my grid is quite rough at the walls. Due to limited computational power, I tried to use as little nodes as possible. I do know that my first cells are around 50 wall distance units large (which is much larger than my turbulence and cfd courses advice me but otherwise the computations will take too long). Maybe this can influence these wall functions so much that it breaks the solver? Although I would expect that the problem then occurs later on in the solver.

FoamFile
{
    format      ascii;
    class       volScalarField;
    location    "0";
    object      alphat;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    Inlet1
    {
        type            calculated;
        value           uniform 0;
    }
    
    Inlet2
    {
        type            calculated;
        value           uniform 0;
    }
    
    Outlet
    {
        type            calculated;
        value           uniform 0;
    }
    
    Wall
    {
        type            compressible::alphatWallFunction;
        Prt             0.85;
        value           uniform 0;
    }
}


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

[dlahaye]

My apologies in case my question that follows comes across as insisting on the previous point raised.

I, however, would like to know what your 0/rho file looks like.

Best wishes, D.

[Jelle Lagerweij]

Oh, I did not set any 0/rho file. Maybe I understood how the package works incorrectly. I thought that setting the thermophysical properties and the initial pressure and temperature correctly would make the code create its own density field from it. I tried to find examples of such initial density state by looking in the tutorials, but couldn't find them. How should such file look like?


(my current 0 folder (for rhoPimpleFoam) holds the following files: alphat, k, nut, omega, p, T, U. It seems that I need to add a rho file, but did I forget others as well?)

[Jelle Lagerweij]

Yes, thanks a lot , I got it now.


Your questions made me think back to the meaning of my inputs enough to find what is wrong.


Before this, I always used incompressible solvers, but now with rhoPimpleFoam and reactingFoam, I switch to compressible solvers. This changes a lot. In my case, I use the thermoPhysical dictionary file to set all the characteristics of my system. This includes making it an ideal gas. Therefore, the density can be computed from the pressure and the temperature.


However, I used gage pressure in my compressible flow simulations. I set the exit pressure to be 0. The ideal gas law will not give a useful result at zero pressure. When the reactingFoam tries to compute some characteristics it will divide by the density and encounter a division by 0...


For users with similar issues. There are at least two major input differences between compressible and incompressible models.


You need absolute pressure instead of gage pressure. Besides that, be carefull, because the incompressible flow is not directly based on pressure. It wants the pressure divided by the density (m^3*Pa/kg). Whereas the compressible solver needs the absolute pressure in Pa.


So if incompressible: use gage pressure/density with dimension of pressure of [0 2 -2 0 0 0 0].


And if it is compressible, set the absolute pressure and give it the as dimension Pa: [1 -1 -2 0 0 0 0].


With kind regards,
Jelle lagerweij


PS: I just finished the rhoPimpleFoam case succesfully and am now implementing my newly earned knowledge in the reactingFoam case. I am really thankful for Domenico's advice as his comments made me think into the right direction.

[dlahaye]

Congrats on the progress that you are making.

Please keep us informed on the progress that you make.

Best wishes, Domenico.

[Jelle Lagerweij]

I ran some small cases last week, while preparing the pre- and post-processing of my larger case. I think that I am nearly done with that and can final start to perform the real model (which I will definitely share over here to show of ). During that progress, however, I encountered some issues which are not major, but annoying.


I use a bash script (se below) to perform all openFoam tasks automatically. The flow is the following, I perform an rhoSimpleFoam model with only air first, before I copy the right files over to the reactingFoam solver, where I start leaking hydrogen. I like to use foamMonitor to check how the solver performs while the solvers are running. However, if I place the line with

foamMonitor -l ...

after the

runApplication $application

line it will only show the convergence plot after the solver is finished. If I have it before, it will break because the postProcessing folder is not created yet. My current workflow is to keep the foamMonitor lines out of my bash file and run them by hand in another terminal during the runtime of the solver. Is there a correct way to automate this? I looked into the ;, && and || options, but I was not that happy with my results.

#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions


# Run the first case
cd SS_initiation
# Set application name
application=$(getApplication)
runApplication surfaceFeatures
runApplication blockMesh

runApplication decomposePar -copyZero
runParallel snappyHexMesh -parallel -overwrite
runApplication reconstructParMesh -constant

rm -r processor*
rm log.decomposePar
runApplication checkMesh

runApplication decomposePar

# foamMonitor -l postProcessing/residuals/0/residuals.dat &
runParallel $application
runApplication reconstructPar
paraFoam
cd ..


# Communicate all the changes and connections between the first and the last case

# Copy the files directly
cp -r SS_initiation/500 Tr_hydrogen/0
cp Tr_hydrogen/0_orig/AIR Tr_hydrogen/0
cp Tr_hydrogen/0_orig/H2 Tr_hydrogen/0
# remove unneeded folder
rm -r Tr_hydrogen/0/uniform

cp -r SS_initiation/constant Tr_hydrogen/constant
rm -r Tr_hydrogen/constant/thermophysicalProperties
cp Tr_hydrogen/constant_old/thermophysicalProperties Tr_hydrogen/constant



# Change time to 0 in 0 folder for ALL files
find Tr_hydrogen/0 -type f -exec sed -i -e 's/location    "500"/location    "0"/' {} \;

# Adjust velocity boundary condition
# start leaking
sed -i 's/uniform (0 0 0);/uniform (0 0 -0.05);/' Tr_hydrogen/0/U


# Run the second case
# Set application name
cd Tr_hydrogen
# foamMonitor -l postProcessing/residuals/0/residuals.dat &
application=$(getApplication)
runApplication decomposePar
runParallel $application
runApplication reconstructPar
reactingFoam -postProcess -latestTime -func yPlus
paraFoam

Additionally, I want to check the mass conservation of hydrogen over my entire control surface. What we need to ascertain with this model is that no accumulation with hydrogen will occur in the room. Therefore, I want to track the leakage (massflow) of hydrogen at my leak inlet, which is easy as it is pure hydrogen and patchIntegrate phi at that inlet will compute that amount perfectly. But I also want to do something similar at the outlet, now I should not integrate the patch over phi, but over x*phi (with x being the ratio of hydrogen). However, I encountered some problems here, phi is a surfaceScalarField and H2 is a volScalarField. How do I interpolate the volScalarField to the outlet patch surface, multiply it with the local flux field and integrate that?


Lastly a short note, I got a weakly turbulent flow and therefore I use kOmegaSST. If I remember it correctly, we would like the wall distance to be between 30 and 300 for that, but my check gives the following result:

yPlus yPlus write:
    writing object yPlus
    patch Wall y+ : min = 0.02835652877, max = 149.8769456, average = 20.64251213

I believe that the wall functions can adjust for the points in the viscous sublayer and buffer layer. Should I track these points with very small y+ values and change my mesh there or is that probably fine? Note that there are some laminar flow regimes as well in my space, so that probably affects these values a lot as well.


Sorry for this long list of questions, but I am getting really close. The last 5 weeks feel like a openFoam crash course to me .

With kind regards,
Jelle Lagerweij

[Jelle Lagerweij]

I got the bash file working as I wanted. I added the following lines before the

runParallel $application

lines:

gnome-terminal --tab --title="monitor2" -- bash -c 'sleep 5 && foamMonitor -l postProcessing/residuals/0/residuals.dat'

Which seems to work just as I want it to.


Now, the only important problem I still have left is this surface integral at the outlet.


With kind regards,
Jelle Lagerweij

[dlahaye]

About the outward flux.

Not sure.

How about multiplying rho * normal-component(U) * species-concentration and integrating this quantity over the outlet patch using appropriate post-processing utilities?