[KingKraut]

Hello alltogether,

I am currently running a simulation with a custom transient solver (based on pisoFoam, version OF2.2.2) on a grid created with cfMesh (v1.1.1, blood vessel geometry with one inlet and several outlets).
The checkMesh output is fine, the option allGeometry gives some errors though. Previous simulations with the same errors and solver did not produce errors, so I assume the mesh is ok.

The simulation runs for a few timesteps, however at some point the adjustable timestepping shrinks until the solver breaks down.
When looking at the results and reconstructing the last timestep I can see a velocity (and pressure, nu, nuSgs, k) divergence right at the inlet of the model, apparently only a few cells on the inlet produce this divergence - and no cells going further in the direction of the flow seem impacted.
Just from looking at it, the cells seem allright (and it is not any of the cells producing the errors in checkMesh -allGeometry - I looked at those with foamToVTK).

I tried running the simulation with different settings, too. But laminar flow (instead of LES), Newtonian transportModel (instead of a custom transportModel) and a ramp to build up the pressure at the inlet linearly (and more slowly) until the oscillating transient pressure is applied did not help.
Also running a simulation with pisoFoam produced the same error at some poinnt.

The only thing that did not produce the error (so far) was applying a velocity boundary condition at the inlet (instead of pressure), however the applied velocities were small in comparison to the resulting velocities with pressure inlet BC. So I don't know if simply the velocity regime is the problem here - if so, I wouldn't know why?

I attached images of the cells on the inlet of the velocity divergence and the checkMesh-file of the grid.

In the following the solver output for the first few timesteps (with pressure ramp on the inlet), this runs fine upt to ~0.9267

Starting time loop

Time = 0.9001

Courant Number mean: 0 max: 0
deltaT = 0.0001
DILUPBiCG: Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Uz, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.009780712546, No Iterations 239
DICPCG: Solving for p, Initial residual = 0.5759537724, Final residual = 0.005631292546, No Iterations 491
DICPCG: Solving for p, Initial residual = 0.4438582217, Final residual = 0.004294747289, No Iterations 492
DICPCG: Solving for p, Initial residual = 0.1389864405, Final residual = 0.001349027544, No Iterations 444
time step continuity errors : sum local = 3.173264302e-08, global = -1.588048215e-09, cumulative = -1.588048215e-09
DICPCG: Solving for p, Initial residual = 0.04888811813, Final residual = 0.0004786009158, No Iterations 446
DICPCG: Solving for p, Initial residual = 0.01723255882, Final residual = 0.0001704062924, No Iterations 479
DICPCG: Solving for p, Initial residual = 0.005970587328, Final residual = 5.818270016e-05, No Iterations 305
DICPCG: Solving for p, Initial residual = 0.002342435853, Final residual = 2.330606533e-05, No Iterations 211
time step continuity errors : sum local = 5.148823953e-10, global = -3.812851673e-11, cumulative = -1.626176732e-09
DICPCG: Solving for p, Initial residual = 0.004126787875, Final residual = 4.07387621e-05, No Iterations 527
DICPCG: Solving for p, Initial residual = 0.004639497045, Final residual = 4.573523558e-05, No Iterations 495
DICPCG: Solving for p, Initial residual = 0.001487301133, Final residual = 1.477091682e-05, No Iterations 482
DICPCG: Solving for p, Initial residual = 0.0006347736756, Final residual = 6.337721064e-06, No Iterations 391
time step continuity errors : sum local = 1.390684203e-10, global = -6.837588022e-12, cumulative = -1.63301432e-09
bounding k, min: -2.384619832e-11 max: 2e-05 average: 7.556430784e-13
inlet pressure: 0.0128
outlet1 flux:2.469759563e-11
outlet2 flux:8.314252151e-11
outlet3 flux:5.459751274e-11
outlet4 flux:5.353705358e-12
outlet5 flux:8.130016969e-12
outlet6 flux:2.426652866e-12
outlet7 flux:5.210336153e-12
outlet8 flux:3.969025875e-12
outlet9 flux:1.351799395e-12
outlet10 flux:1.430630049e-12
outlet11 flux:1.857191531e-11
outlet12 flux:1.722238048e-11
outlet13 flux:1.14786056e-11
outlet14 flux:1.834847085e-11
outlet15 flux:2.094497218e-11
outlet16 flux:1.178583265e-11
outlet17 flux:8.505749459e-12
outlet18 flux:1.513267234e-12
ExecutionTime = 12.49 s ClockTime = 19 s

Time = 0.9002

Courant Number mean: 9.926074585e-06 max: 0.0007871207171
deltaT = 0.0001
DILUPBiCG: Solving for Ux, Initial residual = 0.3379872409, Final residual = 1.940175349e-06, No Iterations 7
DILUPBiCG: Solving for Uy, Initial residual = 0.3505420291, Final residual = 7.908877956e-06, No Iterations 6
DILUPBiCG: Solving for Uz, Initial residual = 0.4578289235, Final residual = 6.114466236e-06, No Iterations 7
DICPCG: Solving for p, Initial residual = 0.3317160549, Final residual = 0.003230101976, No Iterations 244
DICPCG: Solving for p, Initial residual = 0.02041022059, Final residual = 0.0002019589652, No Iterations 616
DICPCG: Solving for p, Initial residual = 0.009536473392, Final residual = 9.372596743e-05, No Iterations 450
DICPCG: Solving for p, Initial residual = 0.002687224167, Final residual = 2.669692757e-05, No Iterations 385
time step continuity errors : sum local = 1.158617378e-09, global = 1.674531243e-11, cumulative = -1.616269008e-09
DICPCG: Solving for p, Initial residual = 0.001564812958, Final residual = 1.523902949e-05, No Iterations 351
DICPCG: Solving for p, Initial residual = 0.002715667817, Final residual = 2.678040409e-05, No Iterations 535
DICPCG: Solving for p, Initial residual = 0.0007602239202, Final residual = 9.778637954e-06, No Iterations 348
DICPCG: Solving for p, Initial residual = 0.0003458572859, Final residual = 9.990685919e-06, No Iterations 27
time step continuity errors : sum local = 4.327538312e-10, global = -4.613208059e-11, cumulative = -1.662401088e-09
DICPCG: Solving for p, Initial residual = 0.0005730344151, Final residual = 9.889904195e-06, No Iterations 484
DICPCG: Solving for p, Initial residual = 0.001180460995, Final residual = 1.178301263e-05, No Iterations 542
DICPCG: Solving for p, Initial residual = 0.0003364532107, Final residual = 9.829179222e-06, No Iterations 247
DICPCG: Solving for p, Initial residual = 0.0001452899322, Final residual = 1.442918761e-06, No Iterations 314
time step continuity errors : sum local = 6.236492828e-11, global = 3.992323544e-12, cumulative = -1.658408765e-09
bounding k, min: -2.150054365e-11 max: 2e-05 average: 6.566971771e-12
inlet pressure: 0.0256
outlet1 flux:7.385637125e-11
outlet2 flux:2.375277029e-10
outlet3 flux:1.579435765e-10
outlet4 flux:1.727871046e-11
outlet5 flux:2.496038986e-11
outlet6 flux:1.464975194e-11
outlet7 flux:1.233799067e-11
outlet8 flux:1.230199031e-11
outlet9 flux:6.136214273e-12
outlet10 flux:3.45799332e-12
outlet11 flux:5.555947686e-11
outlet12 flux:5.217588095e-11
outlet13 flux:3.485109616e-11
outlet14 flux:5.587777564e-11
outlet15 flux:6.342511949e-11
outlet16 flux:3.566514296e-11
outlet17 flux:2.562791476e-11
outlet18 flux:5.371567337e-12
ExecutionTime = 13.9 s ClockTime = 21 s

At some point deltaT starts shrinking, until the simulation breaks down:

Time = 0.927423

Courant Number mean: 0.0008911619681 max: 0.5464201838
deltaT = 7.746316033e-07
DILUPBiCG: Solving for Ux, Initial residual = 0.0003606170857, Final residual = 2.635352112e-06, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.0002949312724, Final residual = 9.535288121e-07, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.000461219898, Final residual = 3.096358811e-06, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.06408172658, Final residual = 0.0006248592931, No Iterations 102
DICPCG: Solving for p, Initial residual = 0.007186542492, Final residual = 7.036256442e-05, No Iterations 465
DICPCG: Solving for p, Initial residual = 0.002839026885, Final residual = 2.767143244e-05, No Iterations 225
DICPCG: Solving for p, Initial residual = 0.0009906055812, Final residual = 9.881980154e-06, No Iterations 298
time step continuity errors : sum local = 1.532491075e-11, global = 3.489345027e-14, cumulative = -2.395069927e-07
DICPCG: Solving for p, Initial residual = 0.00146421456, Final residual = 1.440003767e-05, No Iterations 91
DICPCG: Solving for p, Initial residual = 0.0005728914537, Final residual = 8.980743279e-06, No Iterations 9
DICPCG: Solving for p, Initial residual = 0.000211648722, Final residual = 9.224806136e-06, No Iterations 4
DICPCG: Solving for p, Initial residual = 0.0001106496745, Final residual = 7.280067338e-06, No Iterations 3
time step continuity errors : sum local = 1.132022531e-11, global = 1.697443665e-15, cumulative = -2.39506991e-07
DICPCG: Solving for p, Initial residual = 0.0002004073737, Final residual = 9.546305908e-06, No Iterations 10
DICPCG: Solving for p, Initial residual = 9.019472701e-05, Final residual = 9.983786799e-06, No Iterations 2
DICPCG: Solving for p, Initial residual = 4.343604643e-05, Final residual = 9.865296236e-06, No Iterations 2
DICPCG: Solving for p, Initial residual = 3.00064068e-05, Final residual = 9.86822262e-07, No Iterations 566
time step continuity errors : sum local = 1.534663596e-12, global = -1.36287195e-14, cumulative = -2.395070046e-07
bounding k, min: -21.35773877 max: 160.9365856 average: 0.0004414988092
inlet pressure: 3.510084715
outlet1 flux:5.868566927e-08
outlet2 flux:1.222489716e-07
outlet3 flux:9.922528758e-08
outlet4 flux:1.115984349e-07
outlet5 flux:1.110457612e-07
outlet6 flux:4.707464038e-08
outlet7 flux:1.154937138e-07
outlet8 flux:7.619421753e-08
outlet9 flux:6.283833287e-08
outlet10 flux:7.405393703e-08
outlet11 flux:8.267529338e-08
outlet12 flux:6.300487477e-08
outlet13 flux:1.025960487e-07
outlet14 flux:1.62641933e-07
outlet15 flux:1.692620607e-07
outlet16 flux:7.175510681e-08
outlet17 flux:1.32989048e-07
outlet18 flux:1.475954557e-07
ExecutionTime = 600.19 s ClockTime = 977 s

--> FOAM Warning :
From function Time:perator++()
in file db/Time/Time.C at line 1039
Increased the timePrecision from 6 to 7 to distinguish between timeNames at time 0.9274233115

I am really stuck on this at the moment, I changed the mesh structure (numbers in meshDict), fvSchemes and inputBCs for several variables, but all to no avail. I guess the mesh causes the problems, since all other things worked on different settings and meshs, but I don't understand why this should now work here, since the mesh seems alright (according to checkMesh).

Thanks a lot for anyone looking into this! I hope I did not forget anything of importance! If someone has an idea, what else I could try, please let me know!! All help is highly appreciated!

Thanks again!

Best regards
Johannes

[franjo_j]

Hello,

It is an interesting problem. Can you please try some of the suggestions listed below?

Originally Posted by KingKraut

Hello alltogether,

I am currently running a simulation with a custom transient solver (based on pisoFoam, version OF2.2.2) on a grid created with cfMesh (v1.1.1, blood vessel geometry with one inlet and several outlets).
The checkMesh output is fine, the option allGeometry gives some errors though. Previous simulations with the same errors and solver did not produce errors, so I assume the mesh is ok.

The simulation runs for a few timesteps, however at some point the adjustable timestepping shrinks until the solver breaks down.
When looking at the results and reconstructing the last timestep I can see a velocity (and pressure, nu, nuSgs, k) divergence right at the inlet of the model, apparently only a few cells on the inlet produce this divergence - and no cells going further in the direction of the flow seem impacted.
Just from looking at it, the cells seem allright (and it is not any of the cells producing the errors in checkMesh -allGeometry - I looked at those with foamToVTK).

I tried running the simulation with different settings, too. But laminar flow (instead of LES), Newtonian transportModel (instead of a custom transportModel) and a ramp to build up the pressure at the inlet linearly (and more slowly) until the oscillating transient pressure is applied did not help.
Also running a simulation with pisoFoam produced the same error at some poinnt.

The only thing that did not produce the error (so far) was applying a velocity boundary condition at the inlet (instead of pressure), however the applied velocities were small in comparison to the resulting velocities with pressure inlet BC. So I don't know if simply the velocity regime is the problem here - if so, I wouldn't know why?

Pressure inlet conditions are not the most stable, think of it as an engine pushing the train instead of pulling it. Every obstacle upstream may make the train crash.

In your situation, the obstacles come from the misalignment of the mesh with the flow. It is very likely causes oscillations the of gradient in the neighbour cells that magnify over time.

Can you please make a mesh with two or more boundary layers at the inlet patch? This shall improve mesh-to-flow alignment at the inlet and hopefully solve the problem.

Another suggestion is to activate boundary layer optimisation available in cfMesh.

Can you construct you geometry such that you have a straight stretch at the inlet to develop the flow before entering the domain of interest?

I attached images of the cells on the inlet of the velocity divergence and the checkMesh-file of the grid.

In the following the solver output for the first few timesteps (with pressure ramp on the inlet), this runs fine upt to ~0.9267



At some point deltaT starts shrinking, until the simulation breaks down:



I am really stuck on this at the moment, I changed the mesh structure (numbers in meshDict), fvSchemes and inputBCs for several variables, but all to no avail. I guess the mesh causes the problems, since all other things worked on different settings and meshs, but I don't understand why this should now work here, since the mesh seems alright (according to checkMesh).

Thanks a lot for anyone looking into this! I hope I did not forget anything of importance! If someone has an idea, what else I could try, please let me know!! All help is highly appreciated!

Thanks again!

Best regards
Johannes

It would also be interesting to try various gradient schemes, especially the limited ones, and find out how it behaves then.

Please let me know if any of this helps you solve the problem.

[KingKraut]

Dear franjo_j, dear all,

thanks a lot for your reply, the explanationans and suggestions and please excuse my delayed response.

In fact I did already try the simulations with an extra boundary layer at the inlet and it did not help. The option to optimizeLayers in cfMesh was also switched on for the meshes I created.
(In this regard I asked myself does this also have an effect, if I did not create any boundary layers? In other words does this improve the shape of cells on the boundary in any case? Maybe I will try this sometime...)

I also started some simulation runs with different gradSchemes and divSchemes. It did improve the issue a little bit in that the simulation ran longer. However, it still crashed at some later point for the same reasons...

I attached the fvSchemes from the latest computation. I guess I need to lower the blending factor of the limited divSchemes even more. I chose the schemes according to the mesh quality, and these parameters worked with other hexahedral meshes of comparable quality, where the grid cells were oriented with the flow direction.
Further reducing this factor troubles me a little since this also tampers the solution, which is obtained...

Your idea of changing the goemetry such that I have a straight pipe-like element in the beginning sounds interesting. I will look into that, and hope to be able to do this soon!

Thanks a lot again for looking into this!!

Best regards
Johannes

[KingKraut]

Finally got back to handle this problem now and I am still looking for a solution on this! :-(

By now, I understand, that pressure driven flow is somewhat more complicated to simulate then flow due to a velocity inlet boundary condition. However, application of pressure is the way I need to go, since this is the data I have at hand, and furthermore, I do not want to make any assumption about the velocity profile across the model inlet...

I managed to get the pressure driven simulations to work with a fully hexahedral grid at the inlet oriented with the flow direction (meshed with software ICEM by ANSYS). So I am pretty sure now, that it is the orientation of the cells, that cause the problem. However, meshing with this Ansatz (ICEM) is highly complicated since it requires much user intervention to assure sufficient grid quality. The complexer the geometry gets, the more difficult and time consuming this task becomes - and this forbids its usage on highly complex geometries, which I need to treat also with pressure inlet BC. So fully automated meshing approaches (such as cfMesh or snappyHexMesh) are the thing to do here!
So for these reasons I am dependent on getting these simulations to work.


In the meantime I tried different gradSchemes such as

cellLImited Gauss linear 1;
cellLImited leastSquares 1;

but both only changed the duration until the simulation diverged and broke down.
Currently I am waiting for the performance of the following (I assume most diffusive option) of Gauss linear:

Gauss linear limited 0;

and

faceLimited Gauss linear 1;

But to be honest I am not that optimistic anymore to get things to work with either of this...

So far I did not want to change the divScheme for

div(phi,U) bounded Gauss linearUpwindV grad(U);

to

upwind

because I don't want to make it all even more diffusive, but maybe I have to give this a try?

I attached the latest fvSchemes-dict and the checkMesh-file if anybody is interested to look closer into this. The mesh was created with cfMesh-1.0.1 and according to checkMesh it is totally fine!! the option -allGeometry gives a few bad cells and faces (concaveFaces, lowqualityTetFaces and concaveCells) however all in completely different spots, then where the divergence occurs...

If more information is needed, please let me know. I am glad for anyone looking into this and all advice! As you can see, I am still really stuck in this issue!
Thanks a lot!!



By the way, the software version I am using is OF231 and a custom solver, transportModel and turbulenceModel. However, as I said the errors occur with newtonian transportModel, laminar flow and conventional solver pisoFoam, too.

[KingKraut]

Just tried another option, and extruded the mesh at the inlet for 50 cells (2mm). The first test run again showed the undesired behaviour quite much straightaway...

I will give this possibility another go however with a longer inlet passage of 6mm...

All other ideas are most welcome!!!

[KingKraut]

Furthermore I tried a simulation with the most diffusive divScheme I could think of, upwind. However, again bounding k explodes at some point and time stepping decreases. I guess the simulation will probably break down, soon enough, too...

I am pretty much through with all ideas I have :-(
If anybody else can give another suggestion on how I could solve this problem, I will be very happy to try anything! I don't understand, how the obviously good mesh with checkMesh (better than anything I managed to mesh with ICEM (purely hexahedral and supposedly cells oriented with flow direction) - and there it worked) fails so badly in this case...

Maybe I'll give it a try with snappyHexMesh...

Thanks again for anybody reading this and thinking it through! Please let me know, what you think! All help and comments are highly appreciated!!

[sheaker]

Hello KingKraut.
I have just finish my rans and les pressure driven flow simulations. I had problems with bounding of k and epsilon too! Solution for my problem was setting relaxations factors for k and epsilon to 0.05

In fvSolutions file I have had put:
relaxationFactors
{
U 1;
k 0.05;
epsilon 0.05;
}

to fix my bounding k and epsilon problem.

Hope it helps.

Have a nice day.
Sheaker

[KingKraut]

Hello sheaker,

thank you so much for the response!

I will try this straightaway!

I thought about this option before, however, I was not sure about how "advisable" the usage of underrelaxation is in transient problems?
But my knowledge on this is pretty limited, too...

But thanks again - another straw to clutch! :-)
Have a nice day, too

[KingKraut]

Thanks again for the help, but it did not help either... :-(

[sheaker]

Sorry to hear that.

You said before Your checkmesh looks ok. But I think there are some things You should recheck.

    Mesh non-orthogonality Max: 64.95229193 average: 15.76192355

and

 ***Error in face tets: 33 faces with low quality or negative volume decomposition tets.
  <<Writing 33 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces
    Min/max edge length = 1.947880873e-05 0.0003059821334 OK.
   *There are 738 faces with concave angles between consecutive edges. Max concave angle = 78.5494773 degrees.
  <<Writing 738 faces with concave angles to set concaveFaces

As You mentioned in first post it is probably the mesh problem.

I have just tried simpleFoam on my case and found it diverging due to Uy velocity. There were 1000 iterations of Uy velocity equation so I change relaxation factor for U to 0.5.
In Your case there are hundreds of iteration of pressure equation. For my amateur usage I found pressure most sensitive to non orthogonality and skewness.

Wish You best.

[KingKraut]

Hello sheaker,

thank you very much again for looking deeper into this! This is highly appreciated!

The mesh non-orthogonality max of 64 is quite the best I managed to achieve with different meshing software packages... I am not sure, if I can improve this any further :-/

I looked at the other two checkMesh fails in more detail by means of foamToVTK with paraView and these were spread over the mesh with single cells at very different locations. However, none were on or even close to the inlet, where the discontinuity occurs. So I am not sure if these are the cause... Futhermore, to get rid of these I wouldn't know how to do this.


Another thing that came to my mind over the weekend was the physical and physiologic relevance of my boundary conditions. And it appears to me, that these might actually be quite unrealistic. Possibly the pressure difference between inlet and outlets is too high, resulting in high flow velocites, which thus cause strong gradients too high for this mesh to handle.

However, it still troubles me a lot, that with a mesh of this quality I think I should manage to get a simulation running in some way - at least with the most conservative discretization schemes and relaxation factors.
Like this in the end I should be able to look at the results and say: "yes these boundary conditions are unrealistic, I have to change them."

Still I quite like the idea of playing around with the relaxation factors, and I am not quite done yet adjusting these. Hopefully at some point I will succed to get this running.

Thanks again for the ideas suggestions!

Have a nice day,
Johannes