[cyln]

Hello,

I am doing a steady state incompressible flow simulation with a laminar flow model. My computational domain is a pipe. There are two concentric outlets at the outflow boundary. At the outlets, I have to assign velocity as a massflowrate. I couldnt really find the exact name of this code. So, just to do a test run, i used flowRateInletVelocity with a negative mass flow rate.


1. First of all, what would you recommend for the type? How can I define the mass flow rate at the outlet?

OUTLET1
{
type flowRateInletVelocity;
massFlowRate -10.0;
rho rho;
rhoInlet 1.225;
}

2. When I run the simulation, I get the following error. As I read the threads in the forum, it is caused by the mesh.


#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4 Foam::symGaussSeidelSmoother::smooth(Foam::Field<d ouble>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5 Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#7 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#8 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#9 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 ? in "/opt/openfoam4/platforms/linux64GccDPInt32Opt/bin/simpleFoam"
Floating point exception (core dumped)


3. When I check my mesh using "checkMesh", I get the following error. My mesh is converted from fluent.msh (ICEM) using fluentMeshToFoam. It says Failed 1 Mesh Checks but I cannot really see where the error is. I would be grateful if you could tell me how to interpret this.

Time = 0

Mesh stats
points: 1902
faces: 18200
internal faces: 16704
cells: 8726
faces per cell: 4
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 0
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 8726
polyhedra: 0

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
***Unused points found in the mesh, number unused by faces: 40 number unused by cells: 40
<<Writing 40 unused points to set unusedPoints
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
INLET1 43 32 ok (non-closed singly connected)
INLET2 131 90 ok (non-closed singly connected)
CYLWALL 1143 604 ok (non-closed singly connected)
OUTLET1 42 31 ok (non-closed singly connected)
OUTLET2 137 95 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (0 -0.199986 -0.2) (0.75 0.2 0.199917)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-1.6148e-16 9.11285e-18 2.03217e-17) OK.
Max cell openness = 2.10748e-16 OK.
Max aspect ratio = 4.11226 OK.
Minimum face area = 0.000221586. Maximum face area = 0.00274932. Face area magnitudes OK.
Min volume = 1.83909e-06. Max volume = 4.18896e-05. Total volume = 0.0935029. Cell volumes OK.
Mesh non-orthogonality Max: 51.574 average: 17.1314
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.624006 OK.
Coupled point location match (average 0) OK.

Failed 1 mesh checks.

End


Cheers!

[ChGr]

you boundary condition should work, i guess.

your main problem is the mesh. are you trying to convert a 3D mesh? Try fluent3DMeshToFoam for the conversion...

[cyln]

yea, my mesh is 3D. I used both fluentMeshToFoam and fluent3DMeshToFoam. I get the same error at the end when I check the mesh using checkMesh.

So what is wrong in my mesh?

[ChGr]

not sure with the infos you gave. as a first step i would try to reduce the non orthogonality.

[cyln]

I have reduced the non-orthogonality. Getting the same error unfortunately.

[cyln]

Does anyone else another suggestion?

[alexeym]

Hi,

@cyln

The error, you have posted, can be caused by MANY different things.

Your pressure equation is diverging (since you posted just printStack output, it is just a guess). It can be caused by:

- mesh
- wrong boundary conditions
- discretisation schemes

Since your post provides almost no information (though checkMesh output is surely useful), we need to guess, what is wrong.

1. Did you try PCG instead of smoothSolver?
2. What is you OTHER boundary conditions?
3. What are your discretisation schemes?
4. Are you sure flow is laminar?

[cyln]

Hi alexeym, thank you for your reply,

Yes exactly, my pressure eqn is diverging.

- I tried changing my mesh but I got the same at the end.
- My BCs are as follows:

Velocity:

INLET1
{
type fixedValue;
value uniform (25 0 0);
}
OUTLET1
{
type flowRateInletVelocity;
massFlowRate -3.0;
rho rho;
rhoInlet 1.225;
}
OUTLET2
{
type flowRateInletVelocity;
massFlowRate -6.0;
rho rho;
rhoInlet 1.225;
}
PIPEWALL
{
type slip;
}

-------------------------------------------------------------------------------

Pressure

INLET1
{
type zeroGradient;
}
OUTLET1
{
type fixedValue;
value uniform 0;
}
OUTLET2
{
type fixedValue;
value uniform 0;
}
PIPEWALL
{
type zeroGradient;
}

To me, boundary conditions looks normal. I think flowRateInletVelocity BC at the outlets should also be fine.

- Here are my schemes:

ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss upwind;
div(phi,k) bounded Gauss upwind;
div(phi,epsilon) bounded Gauss upwind;
div(phi,R) bounded Gauss upwind;
div(R) Gauss linear;
div(phi,nuTilda) bounded Gauss upwind;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}

-----------------------------------------------------------------------------------------------------



1. Yes, I have tried PCG and some other solvers instead of smoothSolver. Number of iterations varied bu I got the same error at the end.

4. You are right about the flow. The nature of my flow is turbulent. However, I am running my simulation purposely as a laminar flow. I solved the same problem using CFX and I have to compare my results. OF should still be able to solve it.

[Santiago]

Your system is overdetermined. Put all pressure BC to zeroGradient

Sent from my Lenovo A5000 using CFD Online Forum mobile app

[alexeym]

@cyln

Small additions

1. Rephrasing Santiago in general you cannot arbitrarily fix value of pressure and velocity at a boundary. Since flowRateInletVelocity is a child of fixedValue, you should set pressure BC to something that does not fix value (ex. fixedFluxPressure).

2. Could you also post your SIMPLE dictionary from fvSolution? Do you have non-orthogonal corrector iteration there?

[cyln]

yes, fixedValue pressure was my major mistake. I changed it to zeroGradient. Then it was asking me to define ref cell and pressure, so did I.

- Now it is giving continuity error cuz I can never exactly match the outflow mass flow rate if I do not define mass flow rate at the inflow boundary too. I need to get rid of continuity error while still using fixedValue of velocity at the inlet. I am trying to figure that out.

And here is my SIMPLE dictionary:

SIMPLE
{

pRefCell 0;
pRefValue 0;
nNonOrthogonalCorrectors 10; //0
consistent yes;

residualControl
{
p 1e-2;
U 1e-3;
}
}

[alexeym]

@cyln

1. 10 is too much (though, if you want, you can leave as it is). 2 would be enough. Maybe you could use leastSquares instead of Gauss linear.

2. Since for you case flowRateInletVelocity will do the following:

const scalar avgU = -flowRate_->value(t)/gSum(rho*patch().magSf());

flowRate_, rho, and area of the patch in you case are constants, so you can for example use

postProcess -func 'patchAverage(name=OUTLET1,U)'

to learn exact value of the patch surface and the just set inlet fixed value.

[cyln]

1. I reduced it from 10 to 2.

2. I run the code "postProcess -func 'patchAverage(name=OUTLET1,U)' " to learn exact value of the patch area and set the fixedValue velocity accordingly. Even though I am getting very close value, it is still not exactly the same. See the error below:

Time = 1

smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 1.63795e-05, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 1, Final residual = 0.0268419, No Iterations 1
smoothSolver: Solving for Uz, Initial residual = 1, Final residual = 0.0265898, No Iterations 1


--> FOAM FATAL ERROR:
Continuity error cannot be removed by adjusting the outflow.
Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.
Total flux : 16.6829
Specified mass inflow : 5.84553
Specified mass outflow : 5.8449
Adjustable mass outflow : 1.00136e-18


From function bool Foam::adjustPhi(Foam::surfaceScalarField&, const volVectorField&, Foam::volScalarField&)
in file cfdTools/general/adjustPhi/adjustPhi.C at line 107.

FOAM exiting


Do you suggest me to add the first code in U file?

[Santiago]

Using the so-called 'non orthogonal correctors' does have an impact in your velocity field. Incrementing the intermediate steps on the predictor of PISO/SIMPLE without correcting the velocity generates artificial enstrophy (or "dissipation", as some call it) thus messing the prediction of TKE

Sent from my Lenovo A5000 using CFD Online Forum mobile app

[cyln]

Hi Santiago,

I had the same error even without using that non-orthogonal correctors.
What is your suggestion then?

[Santiago]

Your error is not related to my comment. You can perfectly get a divergence-free velocity field where the TKE is not correct. I just wanted to point out that using the non orthogonal correctors is not inconsequential. Regarding your problem it seems to me that setting the outlets instead of setting the inlet is a bit odd, both from a physical and mathematical perspective. Although numerically is desirable to get an uniform field as an outlet, that's never really the case. Imagine turbulent flow behind a cylinder. Now, on the mathematics part, the equations you are solving are three parabolic and one elliptic, momentum travels unidirectionally, whereas mass conservation is 'omnipresent' this disparity is not well taken by projection methods. Solution of Momentum conservation becomes stiff since information has to travel 'upstream'. Is there any reason you set your exits instead of the inlet?

Sent from my Lenovo A5000 using CFD Online Forum mobile app

[cyln]

I see, thanks a lot for your explanatory reply.

I need to set the inlet velocity BC to 'fixedValue' and outlet velocity BC to 'mass
flow rate'. At the outflow boundary, the velocity is not uniform.

I searched for the name of the BC so that I could set the outlet mass flow rate, but I couldnt find it. So, I did a test run with "flowRateInletVelocity" with a negative mass flow rate.

Anyways, I got stuck with the continuity problem

[Santiago]

You are a FLUENT user, right?

The outflow BC as used in ANSYS does not exist in FOAM, as far as I know. Besides, using such boundaries constraint the physics resolvable by the code. An advective (or even zeroGradient) BC on the fluxes at the outlet is preferred. In the end depends on how much physics are you expecting to see.

Sent from my Lenovo A5000 using CFD Online Forum mobile app

[alexeym]

@cyln

Could you decrease relTol in velocity linear solver? Guess, it is 0.1 now, and this is not enough to satisfy BCs.

Also maybe you can disclose what outflow Bfs did you utilise for CFX simulation?

[cyln]

@alexeym

I varied relTol between 0.1 and 0.01. I am getting the same error.
In CFX, I assigned mass flow rate at the outflow boundary. It is so robust, converge no matter what. I wanna assign the same BC in OF.

@Santiago

I am a CFX user, however, switching to OF. I need to capture the velocity gradient at the outlet as accurate as possible.