[g.freeman]

Hello again,
I am writing here again since OpenFOAM gives me an error and I have no idea how to solve.
I am trying to solve the aerodynamic field around an airfoil with rhoSimpleFoam (I need also the temperature field): the mach number is 0.288.
I prescribed total pressure and temperature at the inlet, static pressure and zero gradient temperature and velocity at the outlet, which seems to be the correct way. By the way, after three time steps, OpenFOAM returns this error:



--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*dFdT)(const scalar) const, scalar (thermo<Thermo, Type>::*limit)(const scalar) const) const
in file /u1/sales/OpenFOAM/OpenFOAM-3.0.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.

FOAM aborting



I tried to change something in fvSchemes file, using limiters, but it doesn't work.

Here you can find my case with the log: naca0012withTemp

Thank you in advance for replying.

[Tobi]

Hello,

the problem here is that OpenFOAM can not calculate the temperature field out of the enthalpy field. If I remember correctly the temperature field is calculated in the thermo - model based on the enthalpy and vice versa. If for any reason the enthalpy / temperature field diverges the result is that message.

For that a) check your simulation before the error; b) check if you underrelax the necessary fields and c) use upwind scheme first.

I can imagine that your mesh causes the problem (some cell).

[g.freeman]

Originally Posted by Tobi

Hello,

the problem here is that OpenFOAM can not calculate the temperature field out of the enthalpy field. If I remember correctly the temperature field is calculated in the thermo - model based on the enthalpy and vice versa. If for any reason the enthalpy / temperature field diverges the result is that message.

For that a) check your simulation before the error; b) check if you underrelax the necessary fields and c) use upwind scheme first.

I can imagine that your mesh causes the problem (some cell).

Hello Tobi,
I did what you suggest and it seems it is all ok; I was using upwind scheme as default. I also changed the mesh (using a coarse and a better one) but the same error appears. With the coarse mesh it appears later, but it still appears.

[Tobi]

Hi,

then there has to be something wrong.
When does that error appear? The best way would be (really) to get the iterations before your solver blow up and check this with paraview. You should be able to see the cells where the problem starts.

[g.freeman]

Originally Posted by Tobi

Hi,

then there has to be something wrong.
When does that error appear? The best way would be (really) to get the iterations before your solver blow up and check this with paraview. You should be able to see the cells where the problem starts.

Hello,
as the log.rhoSimpleFoam reports, the error appears after just 3 iterations:

Time = 3

smoothSolver: Solving for Ux, Initial residual = 0.210427, Final residual = 0.0024659, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.227503, Final residual = 0.00264876, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.990159, Final residual = 0.00557005, No Iterations 1

--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*dFdT)(const scalar) const, scalar (thermo<Thermo, Type>::*limit)(const scalar) const) const
in file /u1/sales/OpenFOAM/OpenFOAM-3.0.x/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) in "/u1/ sales /OpenFOAM/OpenFOAM-3.0.x/platforms/linux64IccDPInt32Opt/lib/libOpenFOAM.so"
#1 Foam::error::abort() in "/u1/ sales /OpenFOAM/OpenFOAM-3.0.x/platforms/linux64IccDPInt32Opt/lib/libOpenFOAM.so"
#2 Foam::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::constTransport<Foam::speci es::thermo<Foam::hConstThermo<Foam:erfectGas<Foa m::specie> >, Foam::sensibleEnthalpy> > > >::calculate() in "/u1/ sales /OpenFOAM/OpenFOAM-3.0.x/platforms/linux64IccDPInt32Opt/lib/libfluidThermophysicalModels.so"
#3 Foam::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::constTransport<Foam::speci es::thermo<Foam::hConstThermo<Foam:erfectGas<Foa m::specie> >, Foam::sensibleEnthalpy> > > >::correct() in "/u1/ sales /OpenFOAM/OpenFOAM-3.0.x/platforms/linux64IccDPInt32Opt/lib/libfluidThermophysicalModels.so"
#4 ? in "/u1/ sales /OpenFOAM/OpenFOAM-3.0.x/platforms/linux64IccDPInt32Opt/bin/rhoSimpleFoam"
#5 __libc_start_main in "/lib64/libc.so.6"
#6 ? in "/u1/ sales /OpenFOAM/OpenFOAM-3.0.x/platforms/linux64IccDPInt32Opt/bin/rhoSimpleFoam"
Abort



I will try to use paraview and undestand which cell creates the problem now.

Thank you for the help

[Tobi]

Hi,

it is just a guess that you have problems on a cell. Maybe it is due to your boundary conditions, then you should find the problem on the boundary face.

Save iteration 1 and 2 and check it out.
Also check your mass conservation and other stuff.

To get a better idea it could help to underrelax with 0.05 to get more iterations (then you will get more data for postprocessing).

[g.freeman]

So, even if I change my underelaxtion factors, it continues to work until t=3, at t=2 the values are in the expected range so it quite hard to detect which cells cause the divergence.
For the boundary condition, I checked which are to use in the book of Ferziger and Peric (computational methods for fluid dynamics)

[Tobi]

Can you share your case? Sometimes is obvious where the problem is located but can not be analyzed without the case (:

Finally, most of the time it is due to cells but in your case it seems that there is something wrong with your BC or solversettings or sth. else.

[g.freeman]

Originally Posted by Tobi

Can you share your case? Sometimes is obvious where the problem is located but can not be analyzed without the case (:

Finally, most of the time it is due to cells but in your case it seems that there is something wrong with your BC or solversettings or sth. else.

Here it is: naca0012withTemp

[Tobi]

Hello,

I think you made a snappyHexMesh meshing and then extrude to 2d?

Your mesh looks okay in the first view but the mesh around your airfoil is very bad (see picture). From one cell to the your layers you have a ratio that is too high.

But that is not the main problem. The problem is related to your boundary conditions. P for inlet and outlet is similar (no pressure drop) and your velocity bc are wrong too. There is no flow and the system is not prescribed in good conditions.

Either, fixed p at outlet and u at inlet or vice versa or set a pressure drop that force the flow...

In addition the mesh should be corrected.

[g.freeman]

Hello again,
first I will refine the mesh around the airfoil and the I'll trying to use the set of boundary condition you wrote.

Thank you for all the help!!

[Tobi]

I do not mean that you have to refine the mesh around your airfoil but the transision between layer and cells are too extreme.

[g.freeman]

Hey all again,
so I changed my mesh in order to correct the AR, and set my boundary condition fixing the pressure at the outlet and the velocity at the inlet.
The solution continues to work for more time step, but soon or later the same error as before appears.
Also the forces are completely wrong respect to that I expect.
Here there is my new case: naca0012heat_new

[Tobi]

The layers are as bad as before.

[opedrofunk]

A few suggestions:
1) In addition to the mesh being a bit funky, your temperature BC's are really weird. Have a look at some of the rhoCentralFoam/rhoSimpleFoam tutorials.
2) try setting nNonOrthogonalCorrectors to around 3.
3) Some of your other BC's are just wrong (using inletOutlet for inlets... etc.; slip conditions for non-vectors...)
4) p BCs seem wrong as well.

Ok, it's pretty clear that your BCs are just not right. Really think through what you are trying to achieve, and try to understand the meaning of each boundary condition.

Also, start from a tutorial case (forwardStep, for example) and use that as a template to start from. Then build complexity incrementally (i.e., use a simple block rather than an airfoil at first and get that right, then add geometric complexity) - I think you'll find that working incrementally toward your goal will be much faster (and less frustrating) than trying lots of stuff and just seeing what sticks...

Good luck!

Kind regards,
Peter

[Tobi]

Are the bc still wrong? Did not check it again

one of the easiest cases to set up

[Carno]

I am getting the same error in rhoSimpleFoam. I tried everything in the discussion.
Please help me. rhoSimpleFoam seems to be extremely picky.

[Tobi]

Dear Carno,

no one knows what you did (you do not provide any information). How should the community help you ?

If you did it like I said above it has to work. Otherwise I think your FOAM is not compiled correct and has errors.

[Carno]

Thanks for the reply.
I am sending information.
It's a MRF problem. A 2m diameter fan is rotating with 250 rad/s. So the tip rotates around 250 m/s.
I used
Inlet BC : total pressure 1 bar
Outlet BC: 1 bar static pressure
Max non-ortho cell is 64 (all tetra)
Kindly guide how to paste code. I can paste the fVschemes and fvSolution here too.
Check Mesh result:
Code:
Mesh stats
points: 612835
faces: 5948384
internal faces: 5438852
cells: 2846809
faces per cell: 4
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 3

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

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
outlet 4494 2406 ok (non-closed singly connected)
walls 332050 166341 ok (non-closed singly connected)
walls_rotor 139415 69893 ok (non-closed singly connected)
walls_hub 29229 14802 ok (non-closed singly connected)
inlet 4344 2331 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (0 -1.02 -1.02) (9 1.02 1.02)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (1.157218e-015 2.861809e-015 -2.727448e-017) OK.
Max cell openness = 2.735671e-016 OK.
Max aspect ratio = 6.275062 OK.
Minimum face area = 3.818967e-006. Maximum face area = 0.005911771. Face area magnitudes OK.
Min volume = 5.255555e-009. Max volume = 0.0001364341. Total volume = 29.37645. Cell volumes OK.
Mesh non-orthogonality Max: 65.20238 average: 16.78288
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.9447552 OK.
Coupled point location match (average 0) OK.

[Carno]

Originally Posted by Tobi

Dear Carno,

no one knows what you did (you do not provide any information). How should the community help you ?

If you did it like I said above it has to work. Otherwise I think your FOAM is not compiled correct and has errors.

FoamFile
{
	version	2.0;
	class	dictionary;
	format	ascii;
	location	"system";
	object	fvSchemes;
}
ddtSchemes
{
	default	steadyState;
}
gradSchemes
{
	default cellLimited leastSquares 1;
	grad(U) cellLimited leastSquares 1;
}
divSchemes
{
	default	bounded Gauss upwind;
	div(phi,U)	bounded Gauss linearUpwind grad(U);
	div(phi,omega)	bounded Gauss upwind;
	div(phi,k)	bounded Gauss upwind;
	div(R)	Gauss linear;
	div((muEff*dev2(T(grad(U)))))	Gauss linear;
}
laplacianSchemes
{
	default Gauss linear limited 0.77;
}
interpolationSchemes
{
	default	linear;
}
snGradSchemes
{
	default limited 0.77;
}
fluxRequired
{
	default	no;
	p	;
}

FoamFile
{
	version	2.0;
	class	dictionary;
	format	ascii;
	location	"system";
	object	fvSolution;
}
solvers
{
	U
	{
		type coupled;
		tolerance	(1.0E-6 1.0E-6 1.0E-6);
		preconditioner	DILU;
		solver	PBiCCCG;
		relTol	(0.1 0.1 0.1);
		maxIter	100;
	}
	omega
	{
		tolerance	1.0E-6;
		preconditioner	DILU;
		solver	PBiCG;
		relTol	0.1;
		//maxIter	100;
	}
	k
	{
		tolerance	1.0E-6;
		preconditioner	DILU;
		solver	PBiCG;
		relTol	0.1;
		//maxIter	100;
	}
	h
	{
		tolerance	1.0E-6;
		preconditioner	DILU;
		solver	PBiCG;
		relTol	0.1;
		//maxIter	100;
	}
	p
	{
		nPostSweeps	1;
		solver	GAMG;
		preconditioner DILU;
		smoother	GaussSeidel;
		nFinalSweeps	0;
		nPreSweeps	0;
		nPostSweeps 2; 
		nCellsInCoarsestLevel	1000;
		cacheAgglomeration	true;
		minIter 3;
		maxIter	100;
		tolerance	1.0E-4;
		agglomerator	faceAreaPair;
		relTol	0.05;
		mergeLevels	1;
	}
}
SIMPLE
{
	nNonOrthogonalCorrectors	3;
	pRefCell	0;
	pRefValue	0.0;
	transonic	true;
	rhoMin	rhoMin [1 -3 0 0 0 0 0] 0.1;
	rhoMax	rhoMax [1 -3 0 0 0 0 0] 10.0;
	residualControl
	{
		U	1.0E-4;
		omega	1.0E-4;
		k	1.0E-4;
		h	1.0E-4;
		p	1.0E-4;
	}
}
relaxationFactors
{
	U	0.05;
	omega	0.05;
	k	0.05;
	rho	0.05;
	h	0.05;
	p	0.05;
}
potentialFlow
{
	nNonOrthogonalCorrectors	10;
	pRefCell	0;
	pRefValue	0.0;
}

FoamFile
{
	version	2.0;
	format	binary;
	class	volScalarField;
	location	"0";
	object	p;
}

dimensions	[1 -1 -2 0 0 0 0];

internalField	uniform 100000.0;

boundaryField
{
	inlet
	{
		type	totalPressure;
		p0	uniform 100000.0;
		gamma	1.4;
		rho	rho;
	}
	outlet
	{
		type	totalPressure;
		p0	uniform 100000.0;
		gamma	1.4;
		rho	rho;
	}
	walls
	{
		type	zeroGradient;
	}
	walls_hub
	{
		type	zeroGradient;
	}
	walls_rotor
	{
		type	zeroGradient;
	}
}