Hi,
I am working on an adapter to couple MBDyn with PreCICE. MBDyn is a moltibody dynamics analysis software. In particular I am working on beam elements. In MBDyn external structural elements are used to exchange load/kinematic information between MBDyn and a CFD solver when beam-like structures are used. Basically I am using the PreCICE interface mesh as a cloud of external points that send forces to the beam nodes and are moved in relation to the beam movement.
I am testing my adapter both with SU2 and OpenFOAM with simple test cases. The interaction with OpenFOAM gives me a lot of problems.
To understand better what is going on, I am working on the vertical flap case with incompressible flow. I made the following steps:

I’ve built a dummy fluid solver which simply reads the interface displacements and writes forces on interface nodes. Using simple concentrated forces (i.e. on the tip) or uniform forces gives the expected shape of the beam. This shows that the basic setup of the MBDyn case and its adapter works.

I’ve also built a dummy structural solver which simply reads the interface mesh the same way of the adapter, reads the forces coming from the flow in the interface (the fluid is simulated with OpenFOAM) and possibly moves the structure in a prescribed way. I have compared the fluid resultants in OpenFOAM and on the adapter and they match. This also tells me that the way the mesh is loaded in the structural solver is correct.

Unfortunately I cannot put everything together. The vertical flap case runs for one time step, converges to a blatantly wrong fluid solution (while the interface mesh barely moves). The fluid solver crashes the next calculation step.

The only thing that I have noticed, between MBDyn and the dummy solver is that, for some unknown reason, without any external forces, the interface mesh shrinks a little bit. I don’t know why, I am investigating with the MBDyn developers.

There is obviously something wrong but I am not able to see it at the moment. all the interface meshes look ok with coherent data.
Sorry for the lengthy post. I could also provide code and or data, but it would have been even lengthier… If anyone is available to exchange ideas and opinions it would be very much appreciated.
Thank you
Claudio

Hi @uekerman,
yes, I am aware, thank you. I started from it as a base. That adapter basically works with plate/shell elements, with no external points mapping: i.e. the nodes of the mbdyn mesh are used as nodes of the interface mesh in PreCICE. I needed to use a slightly different approach.
Claudio

Good steps are typically to test uni-directional coupling

OpenFOAM -> MBDyn, so deformation is not handed back

MBDyn -> OpenFOAM, using some other external force (e.g. a constant force in flow direction) on the structure

For that you have to modify the cplscheme in your config, remove one exchange tag and change to a explicit scheme. See also How to configure one-way coupling

Hi @Makis ,
of course I will. I have just concluded some simulations and I’d like to share some of the results. As I first wanted to validate the adapter, I focused on the Turek-Hron benchmark which is one of the tutorials (this is named FSI3 in the paper). At the moment I am still not able to overcome instability due to mass ratio (solid/fluid) close to 1.
So I decided to consider the other benchmark (named FSI2) which has a mass ratio of 10. My results, considering the tip movement, are the following:

The y displacement has a mean value of 3.35mm with an amplitude 83.9mm and a period of 586ms. The reference has a mean of 1.23mm with an amplitude of 80.6mm.
The x displacement has a mean value of -14.54mm and an amplitude of 11.65mm and a period of 293ms. The reference has a mean -14.56mm and an amplitude of 12.44mm.
I’ve got also other data but I’d like to keep the focus of the post on the configuration aspects and open points:

I need to use a progression of forces: the MBDyn mesh is the same as the PreCICE Solid mesh. The forces on the PreCICE nodes coming form the coupling are passed to MBDyn multiplied by a coefficient:

the first 500ms are used to let the fluid settle, and then for the next 500ms the forces are increased. The slope of the coefficient is really critical for convergence. A steeper slope makes the computation unstable. But this way the simulation requires more time.

I have been using the following coupling parameters:

I can only use Fluid as second participant beacuse otherwise I have issues. But I would be interested in trying other configurations.

I would be interested in discussing some ways to speed-up the convergence and finding a way to overcome the issues regarding the mass ratio.
I am going to update this post when I have new findings
best regards
Claudio

Which internal convergence measures are you using in OpenFOAM and MBDyn? Making those tighter should decrease the number of coupling iterations.

… might be too strict. “Good” values are typically 1e-5 or 1e-6.

The force ramp is a quite standard approach. As an alternative you could apply a ramp to the inflow velocity or re-start from a pre-computed fluid simulation (maybe after 100 timesteps).

Thank you @uekerman , I will try to modify the limit as you suggest.
As to force ramp, I have noticed that in some cases (like this one) the slope can be critical and can lead to divergence even if the simulation would proceed smoothly at “steady state”. I hadn’t much experience with this and it has been a little surprising.
In my case force ramp is much better than inflow velocity ramp, but I wouldn’t say that it is generally true. I am attaching the OpenFOAM configuration files, more or less they are copied from the OF-CalculiX example: fvSchemes.txt (763 Bytes) fvSolution.txt (1.3 KB)

As to MBDyn I am using the following parameters if anyone is familiar with it:

method: ms, .6;
nonlinear solver: newton raphson, modified, 5;
linear solver: naive, colamd, mt, 1;
tolerance: 1e-6;
max iterations: 1000;
derivatives coefficient: 1e-9;
derivatives tolerance: 1e-6;
derivatives max iterations: 100;

Yes, I recommend to use at least a C^1 ramp if you use a second order time integration scheme (e.g. BDF2). For example

0.5 - 0.5 \cos(\frac{\pi}{2} t)

for t \in [0,2]

For “only” a linear ramp, I have also observed problems.

Maybe also interesting here: we are currently working on a proper validation of the FSI3 benchmark with OpenFOAM and dealii (@DavidSCN is in the lead there). Once we are happy with it, we will update the tutorial files.