I have a question about stability in FSI coupling of OpenFoam and CalculiX. I’m currently running a case of a fluid chamber where one side is a deformable membrane. I think one could classify it as an added-mass problem with similar densities.
Question: Why is my case stable with tetrahedron elements for the solid mesh but not with hexahedron elements?
Case description
3D laminar, incompressible (pimpleFoam)
2D scheme of the set up:
For OpenFoam, I use a simple blockmesh. For Calculix, I tried linear Hex elements (He8R) and 2nd order Tet elements (Te10). The meshes have similar numbers of elements.
Graphic of the meshes at the interface: conformal hex-hex mesh. nearest-neighbor mapping used Calculix tet (black), OpenFoam blockmesh (blue). Nearest-projection mapping used
My desired result is the following figure which shows the deformation of the membrane mid-point vs. time. The deformation is similar for both meshes, but the hex-simulation diverges.
What I know so far:
Tet-Elements are stiffer, which leads to the smaller slope in the deformation.
The instability probably lies in the coupling. I’m using implicit methods in the fluid and solid solver (pimpleFoam with Euler, CalculiX dynamic). When simulating the solid part alone, there is no divergence for any element type I used.
The stability of the FSI simulation with hex-mesh is time-step-dependent: The smaller the time step, the earlier the simulation diverges. When using tet-mesh, however, the simulation is stable for a wide range of timesteps (I tested 1e-7 s…1e-4 s).
Lower fluid density helps with convergence.
So… I was ready to experience instabilities at some point, because the setup is tricky, with the fluid only having one inlet/outlet. But I can’t figure out why the tet-mesh configuration is so much more stable than the hex-mesh config.
Ideas:
Element order: The meshes have similar element numbers, but the node number of the tet mesh is higher (due to 2nd order). Maybe this has a positive effect?
Mapping algorithm: For the tet-mesh I’m using nearest-projection, maybe that has an influence on the stability, like some kind of diffusion effect? I would have thought that nearest-neighbor is the “safer” method.
Sorry for the wall of text. I am posting some config and log files below and I can post the full cases upon request! But I hope that I made my question clear for the beginning.
Thanks for the detailed description!
I don’t see any crash or divergence in your files though. Which of the two cases diverges? The one where you have matching meshes? Where do you see the divergence exactly?
Could you please export the coupling meshes in preCICE and upload the vtks. To double-check that you have really matching meshes. Otherwise, the NN mapping can indeed lead to earlier stability problems.
Yes, the case with matching meshes diverges. You’re right, there is no divergence error message. But the CalculiX side gets stuck in calculiation. This is the “never ending simulation” situation that was mentioned in this post. From what I know, this has something to do with divergence, since the residuals and iteration numbers are rising before the “crash” happens.
I attached an archive containing the case files (without logs and solution files). One case with hex-mesh and nearest-neighbor mapping and one case with tet-mesh and nearest-projection mapping.
The vtks do look ok, yes. You can visualize the hex interface mesh using a glyph filter in ParaView.
This is around machine precision. You can interpret the mapping distance as 0. Nothing to worry here. I think this mismatch comes from different treatment of meshes in OpenFOAM and CalculiX. Also the preCICE vtks show a difference in the order of machine precision.
Other than these remarks I am bit out of clues. I could imagine that your case is anyway mathematically ill-defined (and the mapping error for TET meshes allows to hide this a bit better). I don’t understand the OpenFOAM inlet/outlet boundary condition too well. Maybe it helps if you define instead one proper inflow and one proper outflow somewhere? Could be some incompressibility problem. Would at least be interesting to test. You could also try using a compressible solver in OpenFOAM instead.
Otherwise, you could try using a non-matching hex mesh (and NP or RBF mapping in preCICE) and see how the stability behaves then.
Your diagram shows a single inlet/outlet for the water, but it’s unknown whether you are using the OpenFOAM inletOutlet boundary condition there (for velocity). The velocity inletOutlet boundary condition in OpenFoam will set a zero-gradient (Neumann) condition on the velocity if the flow is moving out of your box, and it will set a fixed value (Dirichlet) condition on the velocity if the flow is moving in to your box. You can specify what this fixed value will be. In most cases I’ve seen, it is specified to zero. This results in a boundary condition representing fully-developed flow. Based on your diagram, I don’t think your flow will be fully developed – there doesn’t appear to be enough distance between the membrane and the exit. If it were me, I’d set a fixedValue boundary condition for pressure at that location and use OpenFOAM’s pressureInletOutletVelocity boundary condition for velocity. The pressureInletOutletVelocity boundary condition sets a zero gradient condition for velocity if the flow is moving out of your box (just like the inletOutlet condition), but any part of the boundary has flow coming back into the box (maybe the flow is re-circulating at the edges) the velocity is computed from the fluid flux in the direction normal to the boundary surface.
If it were me, I’d check those BCs and then test out the mapping differences.
Hi @Mike_Tree, thanks for commenting! I should have made my boundary conditions more clear, thanks for pointing that out.
I am in fact using the pressureInletOutletVelocity boundary condition at the inlet/outlet. For the pressure, I chose the totalPressure boundary condition which behaves like fixedValue if the flow is directed out of the domain.
For the sake of completeness, here is the full set of boundary contions for U and p.
p:
internalField uniform 0;
boundaryField
{
interface
{
type zeroGradient;
}
fixedWalls
{
type zeroGradient;
}
bottom
{
type zeroGradient;
}
outlet
{
type totalPressure;
p0 uniform 0;
value uniform 0;
}
}
U:
internalField uniform (0 0 0);
boundaryField
{
interface
{
type movingWallVelocity;
value uniform (0 0 0);
}
fixedWalls
{
type noSlip;
}
bottom
{
type noSlip;
}
outlet
{
type pressureInletOutletVelocity;
value $internalField;
}
}
I think this matches your proposition? If you have further remarks on the boundary conditions let me know!
Edit: P.S.: I also tried using fixedValue with value 0 instead of totalPressure, just in case. It did not change the stability behaviour of this specific case.
Based on this info, I think it’s pretty safe to say that the fluid solver is less likely to be the source of instability. I think an examination of the mapping differences will be very interesting.
I guess more tests are necessary to eventually find the source of this instability.
Since I designed this test case specifically to use the matching mesh for a reference solution, changing the mesh is not what I hoped for.
conformal hex-hex mesh. nearest-neighbor mapping used
Calculix tet (black), OpenFoam blockmesh (blue). Nearest-projection mapping used
