CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

interFoam transport on structured & unstructured meshes

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 29, 2022, 08:58
Default interFoam transport on structured & unstructured meshes
  #1
New Member
 
Join Date: Dec 2021
Posts: 27
Rep Power: 4
finn_amann is on a distinguished road
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<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4 Foam::PBiCGStab::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ??:?
#6 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/finnamann/OpenFOAM/of9release/OpenFOAM-9/platforms/linux64GccDPInt32Opt/bin/interFoam"
#7 Foam::fvMatrix<double>::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
{
    type            zeroGradient;
}
inletSW
{
    type            inletOutlet;
    inletValue      uniform 1;
    value           uniform 1;
}
inletGW
{
    type            inletOutlet;
    inletValue      uniform 1;
    value           uniform 1;
}

U
inletAir
{
    type            zeroGradient;
}
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
{
    type            zeroGradient;
}
inletSW
{
    type            zeroGradient;
}
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;

writeControl    adjustableRunTime;

writeInterval   0.1;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable false;

adjustTimeStep  yes;

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;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
  div(rhoPhi,U)     Gauss linearUpwind grad(U);
  div(phi,alpha)    Gauss vanLeer;
  div(phirb,alpha)  Gauss interfaceCompression;
  div(U) 	        Gauss linear;
  div(alphaPhi.water,C.water) Gauss linearUpwind grad(U);
  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;
    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
File Type: jpg steadystateAlpha.jpg (76.7 KB, 25 views)
File Type: jpg steadystateU.jpg (96.1 KB, 21 views)
File Type: jpg steadystateAlpha_unstr.jpg (84.7 KB, 21 views)
File Type: jpg structuredC.jpg (84.0 KB, 21 views)
File Type: png outline.png (146.5 KB, 17 views)

Last edited by finn_amann; July 1, 2022 at 07:36.
finn_amann is offline   Reply With Quote

Old   July 5, 2022, 11:53
Default
  #2
New Member
 
Join Date: Dec 2021
Posts: 27
Rep Power: 4
finn_amann is on a distinguished road
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<double>(Foam::UList<double> const&) at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/lnInclude/FieldFunctions.C:463 (discriminator 2)
#4 double Foam::gSumSqr<double>(Foam::UList<double> const&, int) at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/lnInclude/FieldFunctions.C:546
#5 Foam::PBiCGStab::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ~/OpenFOAM/OpenFOAM-9/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C:227 (discriminator 1)
#6 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ~/OpenFOAM/OpenFOAM-9/src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.C:170 (discriminator 1)
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/finn/OpenFOAM/OpenFOAM-9/platforms/linux64GccDPInt32Debug/bin/interFoam"
#8 Foam::fvMatrix<double>::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
File Type: jpg dambreak_t0.jpg (196.5 KB, 19 views)
File Type: jpg dambreak_ts53.jpg (197.3 KB, 16 views)
File Type: jpg dambreak_ts54.jpg (197.7 KB, 17 views)
File Type: jpg dambreak_ts107.jpg (195.4 KB, 16 views)

Last edited by finn_amann; July 5, 2022 at 14:19. Reason: clarification
finn_amann is offline   Reply With Quote

Old   July 12, 2022, 12:24
Default
  #3
New Member
 
Join Date: Dec 2021
Posts: 27
Rep Power: 4
finn_amann is on a distinguished road
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
finn_amann is offline   Reply With Quote

Reply

Tags
boundary condition, interfoam, multiphase, transport, unstructured

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Adding diffusion term to interFoam transport equation Gearb0x OpenFOAM Programming & Development 3 February 14, 2023 05:16
Structured and unstructured meshes kveki Main CFD Forum 6 April 5, 2022 00:50
interFOAM: different results for similar meshes bramDeJaegher OpenFOAM Running, Solving & CFD 0 January 31, 2017 19:03
Structured (Cutcell) vs Unstructured (tetrahedral) meshes difference in Fluent result Stan_Mech FLUENT 1 August 20, 2016 21:30
unstructured vs. structured grids Frank Muldoon Main CFD Forum 1 January 5, 1999 11:09


All times are GMT -4. The time now is 17:04.