# interFoam transport on structured & unstructured meshes

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

June 29, 2022, 07:58
interFoam transport on structured & unstructured meshes
#1
New Member

Join Date: Dec 2021
Posts: 27
Rep Power: 4
Hello fellow foamers,

lately, I've been trying out interFoam in conjunction with the phaseScalarTransport functionObject to simulate the transport of a scalar constrained to a single phase (e.g. think of ink in water) that is constantly injected into a steady state through a boundary.

PhaseScalarTransport is part of the foundation version 9 of OpenFoam:
Page in the source code guide
Function object in the github repo

In the following, I'd like to share some findings of how I got things running and which things still don't work. I've attached some images down below that I'm referencing throughout this post.

Generally, my domain is two-dimensional with the following properties:
• Length = 1 m
• Height = 0.3 m
• (Width = 1 m, but it's just one cell so, i. e. the domain is 2D)
The left and right sides of the domain are split up into 3 patches, each being 0.1 m high, to enable different boundary conditions. The "outline" image, that I've attached will give you a little bit more insight.

---

Initially, the lower two thirds of the domain are filled with water (alpha=1) and the upper third is filled with air (alpha=0). The water compartment has an initial flow velocity of 1 m/s in x-direction
(Water flows from left to right by flowRateInletVelocity boundary condition, you can check the boundary conditions for alpha.water, U and p_rgh down below).
From these initial values the simulation ran into a steady state (after about 10 seconds), which can be viewed in the images "steadystateAlpha" and "steadystateU" (and "steadystateAlpha_unstr" to get a glimpse at the unstructured mesh).
Now, this state was reached on a structured mesh (quadratic, ~3000 cells) as well as an unstructured mesh (triangles, ~28000 cells). checkMesh is OK for both grids.

---

Next, I run the simulation anew with the calculated steady state as initial condition and introduce a scalar tracer via the phaseScalarTransport functionObject which is constantly injected over the lower left inflow boundary.
On the structured mesh, this leads to desired results, as the tracer is constantly injected and is constrained to the water phase (see image "structuredC", which is a snapshot at t = 10 s). For me, this already is a big win, which I can definitely build upon in some future studies.

However, on the unstructured mesh, the simulation crashes after about 0.018 s of simulation time with the release-version giving me the classic error message indicating a floating point error (possible division by zero?):

Quote:
 ... Courant Number mean: 0.102808 max: 0.49949 Interface Courant Number mean: 0.00329058 max: 0.36999 deltaT = 0.000352113 Time = 0.0183099 MULES: Solving for alpha.water Phase-1 volume fraction = 0.514418 Min(alpha.water) = -4.29642e-19 Max(alpha.water) = 1.00369 MULES: Solving for alpha.water Phase-1 volume fraction = 0.514418 Min(alpha.water) = -2.90905e-19 Max(alpha.water) = 1.00369 GAMG: Solving for p_rgh, Initial residual = 0.000107181, Final residual = 4.60559e-06, No Iterations 1 GAMG: Solving for p_rgh, Initial residual = 4.67775e-06, Final residual = 1.80101e-07, No Iterations 2 GAMG: Solving for p_rgh, Initial residual = 2.10778e-07, Final residual = 8.19356e-08, No Iterations 1 time step continuity errors : sum local = 4.87583e-08, global = -2.16197e-09, cumulative = 1.44771e-07 GAMG: Solving for p_rgh, Initial residual = 3.32401e-07, Final residual = 5.87953e-08, No Iterations 1 GAMG: Solving for p_rgh, Initial residual = 6.30624e-08, Final residual = 6.30624e-08, No Iterations 0 GAMG: Solving for p_rgh, Initial residual = 6.30624e-08, Final residual = 6.30624e-08, No Iterations 0 time step continuity errors : sum local = 3.75272e-08, global = -2.0845e-09, cumulative = 1.42686e-07 GAMG: Solving for p_rgh, Initial residual = 1.36702e-07, Final residual = 4.08868e-08, No Iterations 1 GAMG: Solving for p_rgh, Initial residual = 4.37115e-08, Final residual = 4.37115e-08, No Iterations 0 GAMGPCG: Solving for p_rgh, Initial residual = 4.37115e-08, Final residual = 4.37115e-08, No Iterations 0 time step continuity errors : sum local = 2.60118e-08, global = -3.05159e-09, cumulative = 1.39634e-07 ExecutionTime = 30.02 s ClockTime = 30 s phaseScalarTransport: Executing phaseScalarTransport: surfaceScalarField alphaPhi.water was not found, so generating it GAMG: Solving for PhiC.water, Initial residual = 0.00212779, Final residual = 4.41967e-05, No Iterations 2 GAMG: Solving for PhiC.water, Initial residual = 5.51592e-05, Final residual = 1.69803e-06, No Iterations 5 GAMGPCG: Solving for PhiC.water, Initial residual = 4.898e-06, Final residual = 1.72034e-08, No Iterations 3 #0 Foam::error::printStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? #2 ? in "/lib/x86_64-linux-gnu/libc.so.6" #3 double Foam::sumProd(Foam::UList const&, Foam::UList const&) at ??:? #4 Foam::PBiCGStab::solve(Foam::Field&, Foam::Field const&, unsigned char) const at ??:? #5 Foam::fvMatrix::solveSegregated(Foam::dict ionary const&) at ??:? #6 Foam::fvMatrix::solve(Foam::dictionary const&) in "/home/finnamann/OpenFOAM/of9release/OpenFOAM-9/platforms/linux64GccDPInt32Opt/bin/interFoam" #7 Foam::fvMatrix::solve(Foam::word const&) at ??:? #8 Foam::functionObjects::phaseScalarTransport::execu te() at ??:? #9 Foam::functionObjects::timeControl::execute() at ??:? #10 Foam::functionObjectList::execute() at ??:? #11 Foam::Time::run() const at Time.C:? #12 ? in "/home/finnamann/OpenFOAM/of9release/OpenFOAM-9/platforms/linux64GccDPInt32Opt/bin/interFoam" #13 ? in "/lib/x86_64-linux-gnu/libc.so.6" #14 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6" #15 ? in "/home/finnamann/OpenFOAM/of9release/OpenFOAM-9/platforms/linux64GccDPInt32Opt/bin/interFoam" Floating point exception (core dumped)
My debug-version gives me another error message entirely, which I don't understand at all, unfortunately. (EDIT: This seems to be connected to a compilation error, as the debug version on my other machine is able to produce an error message similar to the one above. I've opened a GitHub issue concerning this situation https://github.com/OpenFOAM/OpenFOAM-9/issues/12)

Quote:
 ... Starting time loop --> FOAM FATAL IO ERROR: error in IOstream "OSHA1stream.sinkFile_" for operation Ostream& operator<<(Ostream&, const word&) file: OSHA1stream.sinkFile_ at line 0. From function virtual bool Foam::IOstream::check(const char*) const in file db/IOstreams/IOstreams/IOstream.C at line 96. FOAM exiting
I would really like to get the functionObject running on an unstructured mesh. Am I missing something important when it comes to using such a mesh, leading to the simulation crashing?

In the following, I'll give you some more information about my boundary conditions as well as the controlDict, fvSchemes and fvSolution files.

Boundary conditions
The inflow boundary patches (called inletAir (the upper one), inletSW (middle) and inletGW (lower)) on the left have the following boundary conditions:
Code:
```alpha.water
inletAir
{
}
inletSW
{
type            inletOutlet;
inletValue      uniform 1;
value           uniform 1;
}
inletGW
{
type            inletOutlet;
inletValue      uniform 1;
value           uniform 1;
}

U
inletAir
{
}
inletGW
{
type            flowRateInletVelocity;
volumetricFlowRate
{
type            constant;
value           0.1;
}
extrapolateProfile 0;
value           uniform (0 0 0);
}
inletSW
{
type            flowRateInletVelocity;
volumetricFlowRate
{
type            constant;
value           0.1;
}
extrapolateProfile 0;
value           uniform (0 0 0);
}

p_rgh
inletAir
{
type            totalPressure;
p0              uniform 0;
value           uniform 0;
}
inletSW
{
type            fixedFluxPressure;
value           uniform 0;
}
inletGW
{
type            fixedFluxPressure;
value           uniform 0;
}

C.Water (that's the transported scalar)
inletAir
{
}
inletSW
{
}
inletGW
{
type            inletOutlet;
inletValue      uniform 1;
value           uniform 1;
}```
The outflow on the right side is free, meaning zeroGradient for U, alpha.water, C.water and totalPressure, uniform 0 for p_rgh. The lower boundary acts as a wall and the upper as the atmosphere.

controlDict
Code:
```FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
location    "system";
object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     interFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.001;

writeInterval   0.1;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable false;

maxCo           0.5;
maxAlphaCo      0.5;

maxDeltaT       0.01;

functions {
C {
type                phaseScalarTransport;
functionObjectLibs  ("libsolverFunctionObjects.so");
enabled             true;
writeControl        outputTime;
p                   p_rgh;
phase               water; // not needed.
field               C.water;
bounded01           true; // not needed. default is true.
nCorr               0; // set to 2 for unstructured mesh.
D                   0.001;
schemeField         U;
log                 yes;
}
};
// ************************************************************************* //```
fvSchemes
Code:
```FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
location    "system";
object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default Euler;
}

{
default         Gauss linear;
}

divSchemes
{
div(phi,alpha)    Gauss vanLeer;
div(phirb,alpha)  Gauss interfaceCompression;
div(U) 	        Gauss linear;
}

laplacianSchemes
{
default         Gauss linear corrected;
}

interpolationSchemes
{
default         linear;
}

{
default         corrected;
}

fluxRequired
{
default         no;
p_rgh;
pcorr;
alpha1;
}
// ************************************************************************* //```
fvSolution
Code:
```FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
location    "system";
object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{

C.water
{
solver          PBiCGStab;
preconditioner  DILU;
tolerance       1e-08;
relTol          0;
minIter	        1;
}

"pcorr.*"
{
solver          PCG;
preconditioner
{
preconditioner  GAMG;
tolerance       1e-05;
relTol          0;
smoother        DICGaussSeidel;
nPreSweeps      0;
nPostSweeps     2;
nFinestSweeps   2;
cacheAgglomeration false;
nCellsInCoarsestLevel 10;
agglomerator    faceAreaPair;
mergeLevels     1;
}

tolerance       1e-05;
relTol          0;
maxIter         100;
}

p_rgh
{
solver          GAMG;
tolerance       1e-07;
relTol          0.05;
smoother        DIC;
nPreSweeps      0;
nPostSweeps     2;
nFinestSweeps   2;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator    faceAreaPair;
mergeLevels     1;
}

p_rghFinal
{
solver          PCG;
preconditioner
{
preconditioner  GAMG;
tolerance       1e-07;
relTol          0;
nVcycles        2;
smoother        DICGaussSeidel;
nPreSweeps      2;
nPostSweeps     2;
nFinestSweeps   2;
cacheAgglomeration true;
nCellsInCoarsestLevel 10;
agglomerator    faceAreaPair;
mergeLevels     1;
}

tolerance       1e-07;
relTol          0;
maxIter         20;
}

U
{
solver          PBiCG;
preconditioner  DILU;
tolerance       1e-06;
relTol          0;
}

UFinal
{
solver          PBiCG;
preconditioner  DILU;
tolerance       1e-08;
relTol          0;
}

"(k|epsilon)"
{
solver          PBiCG;
preconditioner  DILU;
tolerance       1e-06;
relTol          0;
}

"(k|epsilon)Final"
{
solver          PBiCG;
preconditioner  DILU;
tolerance       1e-08;
relTol          0;
}

alpha.water
{
MULESCorr       no;
nAlphaCorr     	1;
nAlphaSubCycles 2;
cAlpha 		1;
}
}

PIMPLE
{
momentumPredictor	no;
nCorrectors     	3;
nNonOrthogonalCorrectors 2;
pRefCell 		0;
pRefValue 		0;
}```
If there is interest I can also attach my case files, but I have to tidy them up first to make them easily runnable for everyone.

So if someone maybe recognizes some mistakes I might have made, I would be thankful if they could point them out.

Best regards,
Finn
Attached Images
 steadystateAlpha.jpg (76.7 KB, 25 views) steadystateU.jpg (96.1 KB, 21 views) steadystateAlpha_unstr.jpg (84.7 KB, 21 views) structuredC.jpg (84.0 KB, 21 views) outline.png (146.5 KB, 17 views)

Last edited by finn_amann; July 1, 2022 at 06:36.

July 5, 2022, 10:53
#2
New Member

Join Date: Dec 2021
Posts: 27
Rep Power: 4
So I have made a little bit on progress on the analysis of the problem.

I ditched the old case for a dam break scenario, pictured in "dambreak_t0.png". The initial block of water on the left side of the closed domain features a small box of a tracer concentration = 1.

To gain further information, I ran the case with the debug-version of OpenFOAM and wrote out every time step before the crash. Overall, the simulation ran for 110 time steps equaling 0.102005 s before crashing.

However, the values of the concentration start to explode earlier:
From time step 53 to 54, unreasonably high concentration values start to appear at the interface of water and air (I've seen that in my other simulations as well, the interface seems to be the origin of errors each time), see images "dambreak_ts53.png" and "dambreak_ts54.png". Please keep in mind that I didn't adjust the values on the legend in order to keep reasonable visibility.

Another image illustrates the situation at time step 107, where the concetration has spread further and reached values of +-inf, more or less.
---
The console output is also interesting: Until time step 51, my DILUPBiCGStab solver mostly needs 1 to 3 iterations to compute the concentration field, but 52 onwards, it requires around 40 to 50 iterations each time step.

the final error message looks like this:
Quote:
 Courant Number mean: 0.0507017 max: 0.494007 Interface Courant Number mean: 0.00165237 max: 0.494007 deltaT = 0.000752877 Time = 0.102005 MULES: Solving for alpha.water Phase-1 volume fraction = 0.159957 Min(alpha.water) = -1.0303e-18 Max(alpha.water) = 1 MULES: Solving for alpha.water Phase-1 volume fraction = 0.159957 Min(alpha.water) = -1.9697e-18 Max(alpha.water) = 1 GAMG: Solving for p_rgh, Initial residual = 0.0288203, Final residual = 0.000932971, No Iterations 1 GAMG: Solving for p_rgh, Initial residual = 0.00351985, Final residual = 7.765e-05, No Iterations 2 GAMG: Solving for p_rgh, Initial residual = 0.000508791, Final residual = 2.11807e-05, No Iterations 2 time step continuity errors : sum local = 6.58916e-06, global = 2.26358e-07, cumulative = -6.00956e-06 GAMG: Solving for p_rgh, Initial residual = 0.000125085, Final residual = 3.97078e-06, No Iterations 3 GAMG: Solving for p_rgh, Initial residual = 5.75082e-05, Final residual = 2.53905e-06, No Iterations 1 GAMGPCG: Solving for p_rgh, Initial residual = 5.58576e-06, Final residual = 1.99975e-08, No Iterations 3 time step continuity errors : sum local = 6.23871e-09, global = 3.79854e-11, cumulative = -6.00953e-06 ExecutionTime = 186.37 s ClockTime = 187 s phaseScalarTransport: Executing phaseScalarTransport: surfaceScalarField alphaPhi.water was not found, so generating it GAMG: Solving for Phiink.water, Initial residual = 0.159619, Final residual = 0.00244808, No Iterations 2 GAMG: Solving for Phiink.water, Initial residual = 0.0149558, Final residual = 0.000552519, No Iterations 2 GAMGPCG: Solving for Phiink.water, Initial residual = 0.00245973, Final residual = 3.56057e-08, No Iterations 5 DILUPBiCGStab: Solving for ink.water, Initial residual = 2.70413e-05, Final residual = 1.17229e-09, No Iterations 42 #0 Foam::error:rintStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-9/src/OSspecific/POSIX/printStack.C:218 #1 Foam::sigFpe::sigHandler(int) at ~/OpenFOAM/OpenFOAM-9/src/OSspecific/POSIX/signals/sigFpe.C:104 #2 ? in "/lib/x86_64-linux-gnu/libc.so.6" #3 double Foam::sumSqr(Foam::UList const&) at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/lnInclude/FieldFunctions.C:463 (discriminator 2) #4 double Foam::gSumSqr(Foam::UList const&, int) at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/lnInclude/FieldFunctions.C:546 #5 Foam::PBiCGStab::solve(Foam::Field&, Foam::Field const&, unsigned char) const at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C:227 (discriminator 1) #6 Foam::fvMatrix::solveSegregated(Foam::dict ionary const&) at ~/OpenFOAM/OpenFOAM-9/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C:170 (discriminator 1) #7 Foam::fvMatrix::solve(Foam::dictionary const&) in "/home/finn/OpenFOAM/OpenFOAM-9/platforms/linux64GccDPInt32Debug/bin/interFoam" #8 Foam::fvMatrix::solve(Foam::word const&) at ~/OpenFOAM/OpenFOAM-9/src/finiteVolume/lnInclude/fvMatrixSolve.C:326 #9 Foam::functionObjects:haseScalarTransport::execu te() at ~/OpenFOAM/OpenFOAM-9/src/functionObjects/solvers/phaseScalarTransport/phaseScalarTransport.C:434 (discriminator 1) #10 Foam::functionObjects::timeControl::execute() at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/db/functionObjects/timeControl/timeControlFunctionObject.C:119 #11 Foam::functionObjectList::execute() at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.C:684 #12 Foam::Time::run() const at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/db/Time/Time.C:847 #13 Foam:impleControl::run(Foam::Time&) at ~/OpenFOAM/OpenFOAM-9/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl/pimpleControl.C:110 #14 ? in "/home/finn/OpenFOAM/OpenFOAM-9/platforms/linux64GccDPInt32Debug/bin/interFoam" #15 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6" #16 ? in "/home/finn/OpenFOAM/OpenFOAM-9/platforms/linux64GccDPInt32Debug/bin/interFoam" Floating point exception (core dumped)
As mentioned before, it is a floating point error, which seems to stem from a function called sumSqr<double>(Foam::UList<double> const&).

Any ideas why the concentration seems to explode at the interface?

Btw, I created the mesh with GMSH, maybe there's something I have to consider when using gmsh?
Attached Images
 dambreak_t0.jpg (196.5 KB, 18 views) dambreak_ts53.jpg (197.3 KB, 16 views) dambreak_ts54.jpg (197.7 KB, 17 views) dambreak_ts107.jpg (195.4 KB, 15 views)

Last edited by finn_amann; July 5, 2022 at 13:19. Reason: clarification

 July 12, 2022, 11:24 #3 New Member   Join Date: Dec 2021 Posts: 27 Rep Power: 4 Further debugging into the depths of the OpenFOAM code didn't really get me anywhere, since this does not seem to be an error within the implementation or anything like that. But I have since found a thread in which some schemes for scalar transport on meshes with triangular (or tetrahedric? not sure if this is the right word) elements were recommended. I changed the my div schemes to the following: Code: `div(alphaPhi.water,myScalar.water) Gauss limitedLinear01 1;` This way, the values still explode at the interface at some point, but the simulation much better beforehand and more time is simulated. If anyone could input some basic knowledge or experience about simulating transport of a scalar, I would be very grateful . Cheers Finn

 Tags boundary condition, interfoam, multiphase, transport, unstructured

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Gearb0x OpenFOAM Programming & Development 3 February 14, 2023 04:16 kveki Main CFD Forum 6 April 4, 2022 23:50 bramDeJaegher OpenFOAM Running, Solving & CFD 0 January 31, 2017 18:03 Stan_Mech FLUENT 1 August 20, 2016 20:30 Frank Muldoon Main CFD Forum 1 January 5, 1999 10:09

All times are GMT -4. The time now is 09:51.