# twoPhaseEulerFoam and complex geometry

 Register Blogs Members List Search Today's Posts Mark Forums Read

June 17, 2015, 09:40
twoPhaseEulerFoam and complex geometry
#1
New Member

Benoit Soubelet
Join Date: Feb 2015
Posts: 5
Rep Power: 11
Dear all,

I have been struggling recently with twoPhaseEulerFoam in the case of a nuclear fuel bundle with helical wire spacers.
I cannot get the simulation to stabilize and it always ends up with negative absolute temperatures (and void fraction) and then crashes:

Quote:
 PIMPLE: iteration 1 MULES: Solving for alpha.air MULES: Solving for alpha.air alpha.air volume fraction = 4.17513e-06 Min(alpha1) = -1.06435e-16 Max(alpha1) = 0.2 smoothSolver: Solving for e.air, Initial residual = 0.0134295, Final residual = 2.28547e-08, No Iterations 1 smoothSolver: Solving for e.water, Initial residual = 0.0427195, Final residual = 4.8919e-08, No Iterations 1 min T.air 154.616 min T.water 299.883 GAMG: Solving for p, Initial residual = 0.00924996, Final residual = 3.90212e-08, No Iterations 1 GAMG: Solving for p, Initial residual = 8.85205e-05, Final residual = 2.29204e-09, No Iterations 1 PIMPLE: iteration 2 MULES: Solving for alpha.air MULES: Solving for alpha.air alpha.air volume fraction = 4.18606e-06 Min(alpha1) = -1.67499e-16 Max(alpha1) = 0.2 smoothSolver: Solving for e.air, Initial residual = 0.0092352, Final residual = 4.54509e-08, No Iterations 1 smoothSolver: Solving for e.water, Initial residual = 0.0260932, Final residual = 9.17271e-08, No Iterations 1 min T.air 151.052 min T.water 299.87 GAMG: Solving for p, Initial residual = 0.000518398, Final residual = 4.53614e-08, No Iterations 1 GAMG: Solving for p, Initial residual = 6.25783e-06, Final residual = 2.11468e-09, No Iterations 1 PIMPLE: iteration 3 MULES: Solving for alpha.air MULES: Solving for alpha.air alpha.air volume fraction = 4.1868e-06 Min(alpha1) = -1.67499e-16 Max(alpha1) = 0.2 smoothSolver: Solving for e.air, Initial residual = 0.0171473, Final residual = 6.49008e-09, No Iterations 2 smoothSolver: Solving for e.water, Initial residual = 0.000310335, Final residual = 5.80387e-08, No Iterations 1 min T.air -1047.96 min T.water 299.869 GAMG: Solving for p, Initial residual = 0.000622075, Final residual = 6.43898e-08, No Iterations 1 GAMG: Solving for p, Initial residual = 1.0052e-05, Final residual = 5.97372e-09, No Iterations 1 smoothSolver: Solving for epsilonm, Initial residual = 0.0523869, Final residual = 9.31002e-08, No Iterations 1 bounding epsilonm, min: -6091.56 max: 1.47538e+07 average: 32.0509 smoothSolver: Solving for km, Initial residual = 0.720639, Final residual = 1.66572e-09, No Iterations 2 ExecutionTime = 295.76 s Courant Number mean: 9.51017e-07 max: 0.00386195 Max Ur Courant Number = 6.77404 Time = 1.8e-07 PIMPLE: iteration 1 MULES: Solving for alpha.air MULES: Solving for alpha.air alpha.air volume fraction = 4.47712e-06 Min(alpha1) = -1.67497e-16 Max(alpha1) = 0.2 [2] #0 Foam::error:rintStack(Foam::Ostream&) in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so" [2] #1 Foam::sigFpe::sigHandler(int) in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so" [2] #2 [2] at sigaction.c:0 [2] #3 [2] at interp.c:0 [2] #4 Foam:ow(Foam::Field&, Foam::UList const&, double const&) in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so" [2] #5 Foam::tmp > Foam:ow(Foam::GeometricField const&, Foam::dimensioned const&) in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libphaseCompressibleTurbulenceModels.so" [2] #6 Foam::dragModels::SchillerNaumann::CdRe() const in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleEulerianInterfacialModels.so" [2] #7 Foam::dragModel::K() const in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleEulerianInterfacialModels.so" [2] #8 Foam::BlendedInterfacialModel::K( ) const in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleTwoPhaseSystem.so" [2] #9 Foam::twoPhaseSystem::dragCoeff() const in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleTwoPhaseSystem.so" [2] #10 [2] in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/bin/twoPhaseEulerFoam" [2] #11 __libc_start_main in "/lib64/libc.so.6" [2] #12 [2] in "/afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/bin/twoPhaseEulerFoam" [lclrs61:04850] *** Process received signal *** [lclrs61:04850] Signal: Floating point exception (8) [lclrs61:04850] Signal code: (-6) [lclrs61:04850] Failing at address: 0x8c0f000012f2 [lclrs61:04850] [ 0] /lib64/libc.so.6() [0x3cedc326a0] [lclrs61:04850] [ 1] /lib64/libc.so.6(gsignal+0x35) [0x3cedc32625] [lclrs61:04850] [ 2] /lib64/libc.so.6() [0x3cedc326a0] [lclrs61:04850] [ 3] /lib64/libm.so.6() [0x3cee005b38] [lclrs61:04850] [ 4] /afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so(_ZN4Foam3powERNS_5FieldIdEERKNS_5UL istIdEERKd+0x42) [0x7f04806a3ff2] [lclrs61:04850] [ 5] /afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libphaseCompressibleTurbulenceModels.so(_ZN4Foam3p owINS_12fvPatchFieldENS_7volMeshEEENS_3tmpINS_14Ge ometricFieldIdT_T0_EEEERKS7_RKNS_11dimensionedIdEE +0x1bf) [0x7f0483674dbf] [lclrs61:04850] [ 6] /afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleEulerianInterfacialModels.so(_ZNK4F oam10dragModels15SchillerNaumann4CdReEv+0x1c1) [0x7f0482dc74f1] [lclrs61:04850] [ 7] /afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleEulerianInterfacialModels.so(_ZNK4F oam9dragModel1KEv+0x125) [0x7f0482da4315] [lclrs61:04850] [ 8] /afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleTwoPhaseSystem.so(_ZNK4Foam23Blende dInterfacialModelINS_9dragModelEE1KEv+0x4b0) [0x7f0483105900] [lclrs61:04850] [ 9] /afs/psi.ch/project/fast_lrs/workspace/COD/OPFM/OpenFOAM/OpenFOAM-2.3.1/platforms/linux64GccDPOpt/lib/libcompressibleTwoPhaseSystem.so(_ZNK4Foam14twoPha seSystem9dragCoeffEv+0x1e) [0x7f04830f1d1e] [lclrs61:04850] [10] twoPhaseEulerFoam() [0x46b9f2] [lclrs61:04850] [11] /lib64/libc.so.6(__libc_start_main+0xfd) [0x3cedc1ed5d] [lclrs61:04850] [12] twoPhaseEulerFoam() [0x426b91] [lclrs61:04850] *** End of error message *** -------------------------------------------------------------------------- mpirun noticed that process rank 2 with PID 4850 on node lclrs61 exited on signal 8 (Floating point exception). --------------------------------------------------------------------------
Note that the mesh is not perfect, due a very constraining geometry. The checkMesh utility gives me a quite high non-orthogonality (see attached log file).
But I was able to run it with simpleFoam and reach convergence even with turbulence. The results were close to the experimental ones and I thought I could move forward to 2-phase flow.

I worked from the bubbleColumn tutorial, simplifying a little bit the boundary conditions (see attached), but keeping all the other thermo-physical properties.
For the initialization of the turbulent parameters, I use some rules of thumb like for the turbulent intensity L ~ 5% and turbulent length ~ 5% of hydraulic diameter.
So, if I take U=5 m/s as inlet velocity for a given phase.
k = 1/2*(U*I)^2 = 3.1e-2 m2/s2
Then for hydro diam ~ 3e-3m, we have L ~ 1.5e-4m
It gives for epsilon = C_mu^0.75 * k^1.5 / L = 6 m2/s3

I was only able to run a simplified case (without helical wire spacers) in a laminar case and with a "constant" "diameterModel" with "d 1e-4", for the air phase.
But obviously, I'm interested in the full case with turbulence!

I cannot share my mesh, as it is too big. But it is quite straightforward layout: vertical shape, inlet at the bottom, outlet at the top. The side faces are either walls or cyclicAMI boundaries communicating between each other (see picture attached).

In the end, I want to simulate an isothermal flow, so I wish I could get rid of the energy equation. I have the impression it is the source of all my problems.
The several cyclicAMI boundaries might be source of the problem also ... (but again, running fine with simpleFoam).

So if you have ever encountered such issues with OpenFOAM, and found a way around, I would be very happy to learn from your experience!

Best regards
Attached Images
 view.jpg (15.5 KB, 55 views)
Attached Files
 BCs_case.zip (9.0 KB, 10 views) checkMesh.txt (4.1 KB, 7 views) Logfile.txt (9.0 KB, 2 views)

 June 17, 2015, 13:59 #2 Member   Michael Page Join Date: Mar 2009 Location: Quebec, Canada Posts: 36 Rep Power: 17 Hi Benoit I had similar issues with cht simulation. To pass through this error, I first run the simulation in a laminar state for some hundred iterations. After, I restart the simulation with the turbulence activated (you may have to copy your turbulence files in the last timestep) Also, I'm working to make a new solver twoPhaseEulerDyMFoam and I meet some difficulties. Do you have try twoPhaseEulerFoam with cyclicAMI boundary? Best Regards, __________________ Michael Page michael.page@simu-k.com Simu-K inc. www.simu-k.com

 June 17, 2015, 14:29 #3 New Member   Benoit Soubelet Join Date: Feb 2015 Posts: 5 Rep Power: 11 Thanks for the tip I will try this tomorrow morning! To be sure: you are forcing a laminar treatment first, even if your case is turbulent, and then you switch the turbulences on? About you cyclicAMI question, I'm not sure I understand your question. If you ask if my current twoPhaseEulerFoam simulation includes cyclicAMI boundaries, yes it is the case. I have 14 couples of cyclicAMI boundaries.

 June 17, 2015, 14:34 #4 Member   Michael Page Join Date: Mar 2009 Location: Quebec, Canada Posts: 36 Rep Power: 17 Hi, Yes, you set your turbulences as laminar at first. After, switch to your turbulence model. For cycliAMI, you answer to my question Thank you, __________________ Michael Page michael.page@simu-k.com Simu-K inc. www.simu-k.com

 July 2, 2015, 07:40 #5 New Member   Benoit Soubelet Join Date: Feb 2015 Posts: 5 Rep Power: 11 Hello, The turbulence trick didn't work in the end, I was still obtaining crashes after a while. We found a way around by by-passing the solution of these equations in the fvSolution file: Code: ``` "e.*" { solver smoothSolver; smoother symGaussSeidel; //tolerance 1e-7; relTol 2; minIter 0; }``` This way you don't solve at all the energy equations. After several thousands of iterations, I was able to switch them on again with: Code: ``` "e.*" { solver smoothSolver; smoother symGaussSeidel; tolerance 1e-7; relTol 0; minIter 1; }```

 July 2, 2015, 08:17 #6 Member   Michael Page Join Date: Mar 2009 Location: Quebec, Canada Posts: 36 Rep Power: 17 Hi Benoit, You can also play with nAlphaCorr and nAlphaSubCycles. This help in my case. Also, how is your mesh quality with checkMesh? __________________ Michael Page michael.page@simu-k.com Simu-K inc. www.simu-k.com

 July 2, 2015, 08:23 #7 New Member   Benoit Soubelet Join Date: Feb 2015 Posts: 5 Rep Power: 11 Oh I didn't know these options! Thank you, I will look into it. I put the checkMesh output in my first post; according to the utility, the main issue is the non-orthogonality, which is quite high. I'm using a limited 0.33 scheme to partially correct it; I can't go further, it then becomes too unstable.

 July 2, 2015, 10:50 #8 Senior Member   Dongyue Li Join Date: Jun 2012 Location: Beijing, China Posts: 841 Rep Power: 18 Hello, I played with this solver for a long time. Im sure for now its not compatible with tet mesh or moderately bad hex mesh. I tried with 171 - 22x. But from my memory, 23x can not deal with tet mesh either. In all these versions(besides 23x), turbulence model is hard coded. Just check the code. If U are using these versions twoPhaseEulerFoam. I strongly suggest U to use low skewness hex Mesh. Hex mesh is important. If u can not assure your hex's skewness. twoPhaseEulerFoam can only run with laminar. But with tet mesh, it will make delta T extremely small which makes no sense. __________________ My OpenFOAM algorithm website: http://dyfluid.com By far the largest Chinese CFD-based forum: http://www.cfd-china.com/category/6/openfoam We provide lots of clusters to Chinese customers, and we are considering to do business overseas: http://dyfluid.com/DMCmodel.html

 Tags cyclicami, energy, negative temperature, twophaseeulerfoam