|
[Sponsors] |
![]() |
![]() |
#1 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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: Code:
/*---------------------------------------------------------------------------*\ ========= | \\ / 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. Code:
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. Last edited by Jelle Lagerweij; July 22, 2022 at 06:21. |
|
![]() |
![]() |
![]() |
![]() |
#2 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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
HTML Code:
https://openfoamwiki.net/index.php/How_to_add_temperature_to_icoFoam I started to include print statements to see in which line my code broke, it seems to be the following statement where it crashes: Code:
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 |
|
![]() |
![]() |
![]() |
![]() |
#3 |
Senior Member
|
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? |
|
![]() |
![]() |
![]() |
![]() |
#4 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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. Last edited by Jelle Lagerweij; July 26, 2022 at 07:52. Reason: forgot that I used my custom solver. |
|
![]() |
![]() |
![]() |
![]() |
#6 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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: Code:
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 Code:
turbulence->validate(); For N2 it is: Code:
/*--------------------------------*- 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: Code:
/*--------------------------------*- 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.). |
|
![]() |
![]() |
![]() |
![]() |
#7 |
Senior Member
|
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. |
|
![]() |
![]() |
![]() |
![]() |
#8 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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: Code:
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; // } //} // ************************************************************************* // 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. Code:
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; } } // ************************************************************************* // 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. Code:
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; } } // ************************************************************************* // |
|
![]() |
![]() |
![]() |
![]() |
#10 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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?) |
|
![]() |
![]() |
![]() |
![]() |
#11 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
Yes, thanks a lot
![]() 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. |
|
![]() |
![]() |
![]() |
![]() |
#13 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
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
![]() 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 Code:
foamMonitor -l ... Code:
runApplication $application Code:
#!/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: Code:
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 |
|
![]() |
![]() |
![]() |
![]() |
#14 |
New Member
Jelle Lagerweij
Join Date: Jul 2022
Posts: 21
Rep Power: 4 ![]() |
I got the bash file working as I wanted. I added the following lines before the
Code:
runParallel $application Code:
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 |
|
![]() |
![]() |
![]() |
![]() |
#15 |
Senior Member
|
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? |
|
![]() |
![]() |
![]() |
Tags |
ostream error, reactingfoam, sighandler(int) at ??:? |
|
|
![]() |
||||
Thread | Thread Starter | Forum | Replies | Last Post |
Difficulty with U file. Can't run properly <pirun -np 8 simpleFoam -parallel> | cesaranjo | Main CFD Forum | 2 | April 8, 2022 04:00 |