 |
 |
 |
|
|
|
[Sponsors] | ||||
interFoam transport on structured & unstructured meshes |
 |
|
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
June 29, 2022, 08: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:
--- 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:
Quote:
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;
}
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;
}
};
// ************************************************************************* //
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;
}
// ************************************************************************* //
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;
}
So if someone maybe recognizes some mistakes I might have made, I would be thankful if they could point them out. Best regards, Finn Last edited by finn_amann; July 1, 2022 at 07:36. |
|||
|
|
|||
|
July 5, 2022, 11: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:
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? Last edited by finn_amann; July 5, 2022 at 14:19. Reason: clarification |
||
|
|
|
||
|
July 12, 2022, 12: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; 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 |
|
|
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 |
rintStack(Foam::Ostream&) at ~/OpenFOAM/OpenFOAM-9/src/OSspecific/POSIX/printStack.C:218
Linear Mode
