[Saleh Abuhanieh]

Hi Foamers,


I am tying to duplicate the results of a diffuser (Converging Diverging) case, the original results have been obtained using sonicTurboFoam which has been replaced years ago.
In the original case, the author solved the case as steady-state and laminar first with lower outlet pressure then he used the results as initial condition for the transient run.


I tried first to run the case using sonicFoam directly but I receive floating point exception after a number of iterations, then I decided to run the case as steady-state laminar first using rhoSimpleFoam, however, I receive the floating point exception after some iterations but the source of error seems different.


I am using here openFoam 5.0



The two cases are attached without the polyMesh folder due to the attachment size constraint, however, the checkMesh output as follow:


nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 202202
internal points: 0
faces: 400800
internal faces: 198600
cells: 99900
faces per cell: 6
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 99900
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
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
Inlet 90 182 ok (non-closed singly connected)
Outlet 90 182 ok (non-closed singly connected)
UpperWall 1110 2222 ok (non-closed singly connected)
BottomWall 1110 2222 ok (non-closed singly connected)
FrontBack 199800 202202 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-0.17776 0 0) (0.3806 0.066 0.1)
Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
Mesh has 2 solution (non-empty) directions (1 1 0)
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (1.5802e-17 8.90504e-16 5.60963e-15) OK.
Max cell openness = 2.52446e-16 OK.
Max aspect ratio = 5.95444 OK.
Minimum face area = 3.77399e-08. Maximum face area = 0.00019282. Face area magnitudes OK.
Min volume = 3.77399e-09. Max volume = 1.33046e-07. Total volume = 0.00314933. Cell volumes OK.
Mesh non-orthogonality Max: 22.7385 average: 4.27349
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.10047 OK.
Coupled point location match (average 0) OK.

Mesh OK.

End





The mesh seems fine



Any help will be appreciated
Saleh

[Saleh Abuhanieh]

Any ideas?

[allenjohngeorge]

Upload a photo of error you obtained.

[Saleh Abuhanieh]

Hi Allen,


Thanks for your reply.



For the rhoSimpleFoam:


Time = 0.000252

smoothSolver: Solving for Ux, Initial residual = 0.0182727, Final residual = 0.000888658, No Iterations 4
smoothSolver: Solving for Uy, Initial residual = 0.054948, Final residual = 0.00275268, No Iterations 2
smoothSolver: Solving for e, Initial residual = 0.0210301, Final residual = 0.000716406, No Iterations 2
GAMG: Solving for p, Initial residual = 0.260595, Final residual = 0.011883, No Iterations 14
time step continuity errors : sum local = 0.00349148, global = -0.00136865, cumulative = -0.0674022
ExecutionTime = 22.95 s ClockTime = 23 s

Time = 0.000256

smoothSolver: Solving for Ux, Initial residual = 0.0182673, Final residual = 0.00180064, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.0550647, Final residual = 0.00277076, No Iterations 2
smoothSolver: Solving for e, Initial residual = 0.0210038, Final residual = 0.000707412, No Iterations 2
#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::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam:erfectGa s<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ??:?
#4 Foam::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam:erfectGa s<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ??:?
#5 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
#6 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
Floating point exception



For the sonicFoam:


Time = 3.6e-05

Courant Number mean: 0.018126 max: 1.85612
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
smoothSolver: Solving for Ux, Initial residual = 0.0297792, Final residual = 1.78038e-09, No Iterations 1000
smoothSolver: Solving for Uy, Initial residual = 0.189271, Final residual = 9.04104e-08, No Iterations 1000
smoothSolver: Solving for e, Initial residual = 0.104757, Final residual = 7.51345e-05, No Iterations 1000
smoothSolver: Solving for p, Initial residual = 0.0440183, Final residual = 1.27881e-14, No Iterations 1000
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 6.62023e-16, global = 3.55668e-17, cumulative = 9.73302e-16
PIMPLE: iteration 2
smoothSolver: Solving for Ux, Initial residual = 0.0113747, Final residual = 4.76704e-10, No Iterations 1000
smoothSolver: Solving for Uy, Initial residual = 0.110579, Final residual = 9.35679e-09, No Iterations 1000
smoothSolver: Solving for e, Initial residual = 0.0467253, Final residual = 2.49276e-05, No Iterations 1000
smoothSolver: Solving for p, Initial residual = 0.0153356, Final residual = 1.15887e-14, No Iterations 1000
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 6.54113e-16, global = 2.44197e-17, cumulative = 9.97721e-16
smoothSolver: Solving for epsilon, Initial residual = 0.0888693, Final residual = 0.00189198, No Iterations 1000
smoothSolver: Solving for k, Initial residual = 0.0665843, Final residual = 0.000930416, No Iterations 1000
ExecutionTime = 212.68 s ClockTime = 212 s

Time = 4e-05

Courant Number mean: 0.0200068 max: 3.56817
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
PIMPLE: iteration 1
smoothSolver: Solving for Ux, Initial residual = 0.0313989, Final residual = 1.00651e-09, No Iterations 1000
smoothSolver: Solving for Uy, Initial residual = 0.233467, Final residual = 3.06705e-08, No Iterations 1000
smoothSolver: Solving for e, Initial residual = 0.149963, Final residual = 2.32507e-05, No Iterations 1000
#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 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ??:?
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/sonicFoam"
#8 Foam::fvMatrix<double>::solve() in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/sonicFoam"
#9 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/sonicFoam"
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/sonicFoam"



Thank you in advance

[allenjohngeorge]

This happened because your maximum courant number value went above 1 which should be always kept below 1, preferrably below 0.5 . So you have to reduce the timestep ( deltaT ) in controlDict to bring down the courant number value and to get stable solutions.

[Saleh Abuhanieh]

Hi Allen,


Thanks for your reply.



Even if I reduce the time step or make it adjustable, it is just delaying the blowing.

In the original case, the author faced the same problem for the transient, so he was increasing the pressure on the outlet gradually and ran the case as steadystate. So, I am not following the transient run (sonicFoam) any more.



For me, I am not able to get convergence for the steadyState (using rhoSimpleFoam) case.




I think the problem in the BC, I just followed the original case which was produced using sonicTurboFoam.


Any ideas?

Regards,

[peterhess]

Marhaba Saleh,

I was facing such a problem as I used simpleFoam to simulate axial fan in the past and got a solution by adding a reference pressure to the inlet.

Try to add it to the pressure in your case:

Inlet
{
type fixedValue;
value uniform 116800;
refValue $internalField;

}

I am not sure that this suggestion could solve the problem, cause I dont use the solver u r using here, anyway try it and see if you could fix the problem by this suggestion.

Reduce also a relexationFactor for rho field to look like this:

fields
{
p 0.3;
rho 0.01;
}

Give a feed back please, if it works.

I could not test those suggestions on your case, cause there is no blockMesh to generate the mesh.

Regards

Peter

[Saleh Abuhanieh]

Marhaba Peter,


Thanks for your kind reply.


I tried adding the refValue, but that didn't change anything
For reducing the relaxation factor, it is only accelerating the floating point error.


#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::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam:erfectGa s<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ??:?
#4 Foam::hePsiThermo<Foam:siThermo, Foam:ureMixture<Foam::sutherlandTransport<Foam:: species::thermo<Foam::hConstThermo<Foam:erfectGa s<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ??:?
#5 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
#6 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"


Grüße

[peterhess]

OK!

upload please the whole case so I could make some experiments.

Or just add the blockMesh here...

I generated an easy mesh myself and I got a solution without any problem.

That why I am asking for a mesh.

Regards

Peter

[Saleh Abuhanieh]

Hi Peter,


You may find the polyMesh folder in the below link.


https://1drv.ms/f/s!Amq2JoBgpHbgk1-yCmQwwDwnkPMM





If you want to check the original case/results which I am trying to duplicate, you can find it in the thesis of Benjamin Wuethrich the transonic Diffuser case.



Thank you in advance
Saleh

[peterhess]

Moin!

I made some changes in the case and it seams to work.

I will try to summery all the changes I made:

- Adding a reference pressure to the inlet

type fixedValue;
value uniform 116800;
refValue 1e5;


- deltaT 1e-04;

- nNonOrthogonalCorrectors 1;

- relaxationFactors
{
fields
{
p 0.1;
rho 0.01;
}
equations
{
p 0.3;
U 0.9;
e 0.8;
k 0.9;
epsilon 0.9;
}
}


- rhoMin 0.01;
rhoMax 2;


The last one is necessery cause the results shows that the maximum density is 1.5, anyway in your default setups the maximum density is limited by you to be 1.05.

Those changes makes the simulation works.

At least I did not get a divergence yet after 0.3 sec or more than 3000 iterations.

Which change made the simulation works, I am not realy sure, cause I made all of them in one step...

Regards

Peter

[Saleh Abuhanieh]

Hi Peter,


Thank you very much, I got a convergence finally (>20K iterations)

Well, the results not similar to the benchmark one, may be I need to tighten more the residuals control.


If you try the suggestions one at a time, no convergence will occur.
I think it is a mix of relaxation factors and density bounding.
I am sure the refValue doesn't have any effect. I tested it.


Regards,

[Saleh Abuhanieh]

Hi Peter,


Actually, I don't know if we can call this a convergence


When I tighten the residuals, the solver blow out at around 35,000 iterations
One more point, the continuity error keeps increasing


This is the last iteration before the error:


Time = 38167

smoothSolver: Solving for Ux, Initial residual = 0.000524352, Final residual = 7.86366e-06, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.000465268, Final residual = 1.0039e-05, No Iterations 2
smoothSolver: Solving for e, Initial residual = 0.000405433, Final residual = 5.78933e-06, No Iterations 2
GAMG: Solving for p, Initial residual = 2.26562e-05, Final residual = 3.37786e-08, No Iterations 1
GAMG: Solving for p, Initial residual = 2.16707e-05, Final residual = 2.67125e-08, No Iterations 1
time step continuity errors : sum local = 810.117, global = -339.066, cumulative = -1.45143e+07
ExecutionTime = 3084.73 s ClockTime = 3107 s



Regards,

[peterhess]

Hello Saleh!

Well, I did not run the simuilation for more than 3000 iterations. That why I cant see how the progress is going forwared...

Anyway, The results in paraView für the 3000 iterations where actualy plausible and the reiduals where realy fine, as I uploaded im my last comment...

I suggest to compare the results between 3000 and 30000 iterations using paraView!

The 3000 iterations are realy somth and nice, that why I suggest to use it as a compare stage...

If the results are the same or almost the same, then the results are acceptable...

If not, then something is going wrong...

That the residuals sometimes fluctuating so much is not nice but acceptable...

Addetionaly, the Final residual are realy low...

I would accept the results, if the 3000 and 30000 iterations results are identical...

Regards

Peter

[Saleh Abuhanieh]

Hi Peter,


Actually there is difference between the results from 3000 to 30000 iterations and both don't match the reference/benchmark results.


You may check the attached screenshots.


What about the high continuity error?



Regards,

[peterhess]

Hello Saleh,

The main question ist, are the results of 3000 and 30000 identical?

If yes, we can say, that the solution is under its conditions acceptable, for itself

After that I suggest trying to do the next step, comparing the setups of the original simulation with yours...

As example, try to be sure that the pressure at the inlet and the outlet are realy statical as you defind them and not totalPressure as I gess...

By the way, I would also put a question mark on the results of the original simulation. They are not Quran

Regards

Peter

PS: I see they are not identical... I will run the simulation myself and have a look

[Saleh Abuhanieh]

Hello Peter,


As you noted .. they don't match.
All the pressure values are static and the case setup is identical, however those results have been produced using sonicTurboFoam, so may be I need not to follow the same setup .


The results in the original case have been validated using experimental data and by using another CFD code, so I trust his result (more than mine at least )



Thank you again for your support

[peterhess]

Originally Posted by Saleh Abuhanieh

Hello Peter,

All the pressure values are static and the case setup is identical, however those results have been produced using sonicTurboFoam, so may be I need not to follow the same setup .

Well, I compare the results using rhoSimpleFoam...

As I mentioned before, I did not use sonicTurboFoam in the past and I cant comment about it...

All suggestions made above are for rhoSimpleFoam.

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

Three more notices...

- The bottom wall is defined to be wall!

Is that right or you mean symmetrical?

If you want to simulate a nozzle and you are simulating the half, then it must be type symmetrical!

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

- Re number?

in your case:

rho at the inlet = 1.5 (I gess)

charasterestical length is the height of the tunnel (nozzle) X 2 (by defination in heat atlas) = 2 X 0.066 (supposing the bottom wall is wall)

U = 150.6

Mu = 1.8 e-5

Re = 1.5 X 2 X 0.066 X 150.6 / 1.8e-5

Re = 33 e6

I would say the flow is turbulent (I gess) and you are simulating laminar!

And one more time I am talking about rhoSimpleFoam here...

Is the original simulation laminar or turbulent?

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

If you have experiment results, are the front and back, as you defined in the case realy empty or they should be walls?



I still simulating. I will give a feed back later.

Regards

Peter

[peterhess]

Hello Saleh!

I executed the simulation with rhoSimpleFoam applying the changes I suggested above and used OF6.

I got stable residuals until about 14000 iterations, then the simulation became unstable, until it explode by 45000 iterations.


The divergence happens very late cause the relaxationsFactor of the pressure has been reduced. The simulation became very conservative like that...


Regards

Peter

[Saleh Abuhanieh]

Originally Posted by peterhess

Well, I compare the results using rhoSimpleFoam...

As I mentioned before, I did not use sonicTurboFoam in the past and I cant comment about it...

All suggestions made above are for rhoSimpleFoam.

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

Three more notices...

- The bottom wall is defined to be wall!

Is that right or you mean symmetrical?

If you want to simulate a nozzle and you are simulating the half, then it must be type symmetrical!

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

- Re number?

in your case:

rho at the inlet = 1.5 (I gess)

charasterestical length is the height of the tunnel (nozzle) X 2 (by defination in heat atlas) = 2 X 0.066 (supposing the bottom wall is wall)

U = 150.6

Mu = 1.8 e-5

Re = 1.5 X 2 X 0.066 X 150.6 / 1.8e-5

Re = 33 e6

I would say the flow is turbulent (I gess) and you are simulating laminar!

And one more time I am talking about rhoSimpleFoam here...

Is the original simulation laminar or turbulent?

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

If you have experiment results, are the front and back, as you defined in the case realy empty or they should be walls?



I still simulating. I will give a feed back later.

Regards

Peter

Hi Peter,


- the case is not symmetrical, it is a convergent/divergent channel with a flat bottom


- As I mentioned in my first post, the case is turbulent and transient, however, we run the case as steadyState and laminar with lower outlet pressure to initialize the transient case (I am here now). I didn't get convergence directly (without initialization from steadyState) neither the author of the original case.


- The case is pure 2D, the third direction is empty


You may check the case in section 4.2 (page 65/100)


https://www.research-collection.ethz.ch/bitstream/handle/20.500.11850/150401/eth-30357-01.pdf