Divergent FSI simulation with Reynold's number = 100 and above

Hi,
I have been trying to reproduce TurekHron FSI benchmarking results.
It involves a parabolic inlet condition.
And FSI simulation works fine for Re = 20, but gives the following error for Re=100.


I am using OpenFOAM and CalculiX for the FSI simulation coupled with preCICE.

I have attached my log files for both the successful and unsuccessful cases here. Both Cases

I am changing my kinematic viscosity to change Re.
I am stuck and any help would be really appreciated

What you experience is a divergence:

PIMPLE: Iteration 11
smoothSolver:  Solving for Ux, Initial residual = 0.7320385686, Final residual = 20.08690705, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.7140772146, Final residual = 5.819994424, No Iterations 1000

Look at your final residual: it is too high (it should be much lower than 1). A hint for that is that the simulation has reached the maximum number of iterations (set here to 1000). This means that the results you get are very unphysical.

preCICE also has troubles converging at this situation, already from the first time step (notice the it 50 of 50 part):

PIMPLE: Iteration 50
smoothSolver:  Solving for Ux, Initial residual = 6.73806608e-16, Final residual = 6.73806608e-16, No Iterations 0
smoothSolver:  Solving for Uy, Initial residual = 2.411149517e-17, Final residual = 2.411149517e-17, No Iterations 0
GAMG:  Solving for p, Initial residual = 4.206542465e-11, Final residual = 4.206542465e-11, No Iterations 0
time step continuity errors : sum local = 7.039356843e-16, global = 6.156768306e-18, cumulative = -5.143514102e-12
GAMG:  Solving for p, Initial residual = 4.206542214e-11, Final residual = 4.206542214e-11, No Iterations 0
time step continuity errors : sum local = 7.039356843e-16, global = 6.156768306e-18, cumulative = -5.143507945e-12
ExecutionTime = 88.09 s  ClockTime = 92 s

---[preciceAdapter] [DEBUG] Writing coupling data...
---[preciceAdapter] [DEBUG] Advancing preCICE...
(0) 03:39:37 [mapping::PetRadialBasisFctMapping]:873 in printMappingInfo: Mapping Forces0 conservative from Fluid-Mesh-Faces (ID 0) to Calculix_Mesh (ID 2) for dimension 0) with polynomial set to separate
(0) 03:39:37 [mapping::PetRadialBasisFctMapping]:873 in printMappingInfo: Mapping Forces0 conservative from Fluid-Mesh-Faces (ID 0) to Calculix_Mesh (ID 2) for dimension 1) with polynomial set to separate
(0) 03:39:37 [mapping::PetRadialBasisFctMapping]:873 in printMappingInfo: Mapping Forces0 conservative from Fluid-Mesh-Faces (ID 0) to Calculix_Mesh (ID 2) for dimension 2) with polynomial set to separate
(0) 03:39:37 [mapping::PetRadialBasisFctMapping]:873 in printMappingInfo: Mapping Displacements0 consistent from Calculix_Mesh (ID 2) to Fluid-Mesh-Nodes (ID 1) for dimension 0) with polynomial set to separate
(0) 03:39:37 [mapping::PetRadialBasisFctMapping]:873 in printMappingInfo: Mapping Displacements0 consistent from Calculix_Mesh (ID 2) to Fluid-Mesh-Nodes (ID 1) for dimension 1) with polynomial set to separate
(0) 03:39:37 [mapping::PetRadialBasisFctMapping]:873 in printMappingInfo: Mapping Displacements0 consistent from Calculix_Mesh (ID 2) to Fluid-Mesh-Nodes (ID 1) for dimension 2) with polynomial set to separate
(0) 03:39:37 [impl::SolverInterfaceImpl]:390 in advance: it 50 of 50 | dt# 1 | t 0 of 10 | dt 0.001 | max dt 0.001 | ongoing yes | dt complete no | read-iteration-checkpoint | 

Furthermore, if we look at the succeeded case with Re=20, something looks very odd already in the inlet boundary (so, a CFD-specific problem, nothing with coupling yet):

I would recommend that you take a step back and run the CFD model alone for an extended simulation time and see that everything is fine there first. Eitherway, the last result files I see in the succeeded case are for t=0.6s. What happens if you run it for longer?


Parabolic inlet profile

I am not sure what is going wrong here:

inlet
    {
        type fixedProfile;
        profile csvFile;
        profileCoeffs
        {
          nHeaderLine 0;
          refColumn 1;
          componentColumns 3(2 3 4);
          separator ",";
          mergeSeparators 0;
          file "Fluid/0.orig/0.2U.csv";
        }
        direction (0 1 0);
        origin 0,0,0;
    }

but it would be a good idea to try with groovyBC (part of swak4Foam). I don’t have a config with me at the moment, but if you make one (or find the problem with the csv approach), it would be nice to post it here.

Keep in mind that you also need a time-varying velocity inlet.


PETSc error

PETSc just catches and forwards the errors/signals of any application using it (same with MPI). This indeed makes it very cryptic. The error here is: Caught signal number 8 FPE: Floating Point Exception,probably divide by zero and it usually means a numerical explosion somewhere (in this case, in the velocity calculations).

Hey,
I have corrected the parabolic input error.
For some reason I am unable to get groovyBC to run. So, I opted for a CSV file input.

And I think I have corrected my case taking your remarks.

This is my test case for the 3rd FSI benchmarking of Turek Hron
Only changes are
1)I am taking U mean = 1 msec-1 instaed of 2 msec-1
2) I am correcting the Re to 200 by changing viscocity
3) I am accounting for the Ae given in the paper by changing E* (Young’s Modulus)
4) Without giving a gradual velocity input, I am giving a developed flow from the beginning

I am enclosing the case here with the log files. FSI3

The residuals are way below 1 this time. Still the simulation halts.
I do not understand the error at all.

Thank you in advance.
P.S I realize it’s a Sunday and I am really grateful for your help.

Quite some progress, great! :smiley:

Now the problem is in the Solid participant, as I see in the Solid.log:

*ERROR: solution seems to diverge; please try 
 automatic incrementation; program stops
 best solution and residuals are in the frd file

I cannot look into more details at the moment, but maybe you get a trail. I would open the CalculiX results in CGX and see how they evolve. Is your structure maybe too light (so you would need to adjust the coupling parameters)?

You are already using implicit coupling and RBF mapping, so this should be fine.

The error you see on the Fluid side just means that the other participant has crashed and the Fluid does not find any data to receive:

Receive failed: read: End of file

Hello.
The CalculiX part was causing the problem.
I think I solved it for the most part.
On top of that, I used a groovyBC instead of a csv file to give the parabolic inlet flow.
The steady flow test case FSI1 fully ran and gave perfect results when compared to the benchmark. It can be found here.

While, the unsteady state for Re = 100 was properly set and giving nice results until the time step reached 3.589.
That’s when pimpleFOAM started crashing. The error is in the picture.

The log files, and Fluid.foam and whole case is attached here.
I think I am soo close. It was perfect until the point it ran.
I am attaching both the unsteady case here.
Any help to further develop the case would be much appreciated.

Thank you.

You get a max Courant number of 5. It could be worthwhile trying a smaller timestep size.

I tried that but the simulation is taking too long.
So, then I tried using adjustable time step, making use of subcycling of precice when required.
Same result.

Although, by changing the solver from GAMG to PCG… The simulation was stable for a longer time, 8 seconds this time.
Do you have any idea what it means or how it can be improved?

So, this sounds like another indicator that the problem is the fluid convergence. I am no OpenFOAM expert, but GAMG sounds like a algebraic multigrid solver and PCG like a conjugate gradient solver. But the fact that exchanging both has a crucial impact on the robustness shows that you are already at the stability limit.

If making the timestep size smaller is too expensive you could try to make the fluid mesh coarser (will have the same effect on the Courant number).

Also, other than exchanging the solver you could also try to play around with the solver convergence limits.

Hello @MakisH, @uekerman

I finally figured it out.


It’s stable. The y-displacement plot of the tip is as shown which properly aligns with the paper with a minor error.
Thanks for all the help.

It was difficult to replicate the paper as a lot of forum discussions outside are a little misguiding. Many modifications for the solid as well as fluid were required
Anyone who is newly working with preCICE and see the modified Turek-Hron case in the tutorials repository and wants to replicate the original case might be facing a similar difficulty.
I don’t want to get ahead of myself but is it possible that I send the case after refining it a little more and it can serve as a standard tutorial on the preCICE github tutorial repository?

Thanks,
Nithin.

That’s great, @nithinadidela

We would be very happy if you could open a PR in the tutorials repository (please let us know if you need help with this). We are currently trying to validate the Turek case ourselves as well and wanting to update the case.

Hey @uekerman

Sorry for the delay.
I finally found time and opened a PR in the tutorials repository. I think I did it right. I am new to this. So, it’d be great if you can inform me if I messed up something.

Thank you.

Thank you very much, @nithinadidela! I don’t see any PR from you in precice/tutorials, could you please post a link here?

Hi @MakisH
This might sound funny but this is the first time i am using github

Here’s the link PR

It is a completely normal situation for our users to not be familiar with GitHub. We are here to bridge computer science with mechanical engineering! :wink:

What you currently have is a branch in your fork, so you are ready to open a pull request! You can follow this guide to open a Pull Request to the original repository. We will then try to run it, go through the code, and probably give you some feedback to improve the contribution before merging it. This is what we also do among us in our daily work.

Hints for a good PR:

  • Write a description on what this does, any known limitations/problems, and maybe comment on special aspects of the code
  • Give some starting point for the review: is there anything you would like us to check?

It may sound a bit complicated, but it is actually a very simple process through which one can learn a lot!

Hi @MakisH
I think I did it right this time.
Please let me know if I messed up again.

Here is the link for the pull request 64

Perfect! Thank you! :smiley:
Let’s continue the discussion there.

Hi @MakisH
Are the cases running?

Also I do have one more case Lid Driven Cavity with time varying cosine input. It’s validated with a paper. Should I send a PR with that case? Would it be helpful for preCICE users?

Hi @nithinadidela,

I will have a look at your cases hopefully next week. We are currently a bit swamped with the new preCICE release.
The driven cavity case sounds also very interesting. Yes, please also open a separate PR for this case. Thanks a lot, appreciated!