[Zymon]

Hi,

I need your help. I am dealing with this problem for a couple of weeks.

I am simulating a membrane pump. A pressure plate presses a membrane into the fluid chamber. Due to the reduction in volume the blood is pressed out of the chamber. The membrane has a low Young's Modulus (about 5 MPa) and its density is almost the same as the one of blood. So the Membrane is very flexible, and here is the main problem of the FSI simulation, I guess.

With a hexahedral SOLID185 mesh, frictionless Augmented Lagrange Contact Conditions (FKN: 0.1, agressive Update each iteration, Auto Detection Value of Pinball Region) and a direct solver the uncoupled FEM simulation (ANSYS transient structural, dummy pressure at the FS-interface) converges well. Coupling it with CFX it diverges however in the first time step. Specifically in the FEM part of the second coupling iteration. I added a plot of the convergence history.

So, when I use a much higher Young's Modulus and density of the membrane the FSI simulation converges (untill the point the simulation crashs because of unphysical material behaviour). Using a tetrahedral FEM mesh (also SOLID185), pure penalty algorithm and lower FKN value the FSI simulation runs for several time steps even with not so high values (the standalone FEM case has problems after some time steps with this configuration). However, using a lower time step the simulation diverges again.

I expected an artificial added mass instability because of the density ratio and the general behaviour. So I tried several things:
- very low underrelaxation factors
- source term at the interface with a very high mass flux pressure coefficient (up to 1e37). This lowers the force at the interface (can be seen in a monitor) but has no effect on convergence
- double precision
- higher time steps
- and so on. All the tipps I could find in this forum.

But it does not help. The convergence history always looks almost the same and the simulation diverges in the first time step. Has anybody experience the same problem? Are there any other parameters I could apply?

[ghorrocks]

I am just guessing but your light and flexible membrane will be difficult to converge as it will oscillate with flow. Have you modelled this as a dynamic model in the FEA? Have you included the membrane mass? You might need to increase the membrane mass to dampen these oscillations.

[Zymon]

Now I found out that the CNCHECK,AUTO Command in transient structural is killing any FSI simulation I tested. However, for the FEM simulation it is quite helpful. So, switching to CNCHEK,ADJUST the FSI simulation runs for some more coupling iterations untill it crash. But now it ends with the error "One or more elements have become highly distorted" (see first picture). Looking at the force acting on the interface a sudden peak can be seen (see second picture). This indicates an artificial added mass instability. So I try to change the value of the source term and the underrelaxation factors again.

Thank you ghorrocks for your reply. I thought about oscillations of the membrane as well. This can also be the reason for the distortion of the elements. How can I add mass to the membrane? I am importing the geometry into the Design Modeler and adding material properties with the engineering data. Should I include any damping factor (beta, what value?)?

[Zymon]

Ok, short update. I increased the mass flux pressure coefficient to 1e5 and also set the relaxation factors to reasonable values of 0.6 (higher value as before). Now the force at the interface stays below 1 during the first time step. The second time step looks quite good right now as well. So the source term has an impact now. The CNCHECK,AUTO command killed this effect, so it should not be used when perfoming FSI simulations!

Maybe this helps someone. I keep you updated whether the simulation runs to the end or not.

[ghorrocks]

You should set the mass and damping to match your material properties.

Alternately, it is likely that the mass and/or damping is so small as to be irrelevant, but some mass and/or damping is required for numerical stability. In this case do a sensitivity analysis to sweep through a range and find a value which is high enough to be stable but small enough to not affect the results.

[Zymon]

Where can I set the mass? And should I use the k-Matrix Damping Multiplier for the Membrane? Thanks a lot!

[stumpy]

The mass flux pressure coefficient of 1e5 sounds way too high from my experience. It's possible this is correct, since the size of this coefficient depends on the size of the other coefficients in the matrix, which you really don't know. Could you post your CFX CCL from the .out file?

[Zymon]

Thanks a lot stumpy for your help stumpy. I have seen from the other posts that you have a good knowledge of these kinds of problems. I attached the ccl file to this post.

A mass flux pressure coefficient of 1e5 is indeed too high, I guess. The simulation converges quite well, but the first results don't look promising. I get increased velocities of the flow at the interface that are very unphysical. In the rest of the domain the flow is almost stationary. Using a value of 10e-3 or lower, the simulation gets unstable again. Furthermore, now the RMS P-Mass residual does not converge in each CFD loop. What do you think about adding a mass and damping factor?

[Zymon]

What about using a lower mass flux pressure coefficient in combination with a lower relaxation factor for the force? Maybe there is an optimal compromise between stability, accuracy and computational effort.

[stumpy]

A few points to check...
- The Interface is the only moving boundary, which means it's adjacent to stationary boundaries. This means in Mechanical the edges of your FSI interface should be adjacent to fixed supports or constrained in some way to stop them moving. I just mention this since it's a common setup error.
- I assume the structure is not pre-stressed. Your initial conditions are 1[mm Hg] for pressure. Just be aware that you are simulating the case of an unstressed structure that suddenly (in 0.001[s]) sees a step pressure change of 1[mm Hg]. That might be OK, since 1[mm Hg] is not too big, but it would be better to use a reference pressure of 1[mm Hg] and and opening pressure of 0 to avoid that bump at the start of the transient.
- I would set the (explicit) Under Relaxation Factor to 1. The mass flux pressure coefficient is a form of (implicit) under-relaxation, so it's easier to just have 1 setting to change, and it will converge faster without using explicit under relaxation
- You certainly shouldn't need 40 coupling iterations, but hopefully the convergence criteria will take care of using fewer coupling iteration.
- Use a Maximum Number of Coefficient Loops of 4. One way to think of the mass flux pressure coefficient is to think of it as slowing down the pressure response at the interface. If you perform enough coefficient loops it's still going to end up in the same place, it just takes longer to get there. So the balance should be between the number of coefficient loops and the coefficient value. A large coefficient and lots of loops gets you to the same place as a smaller coefficient and fewer loop. You may as well fix the number of loops to a small value then tune your coefficient.
- After making the above changes set your coefficient to a small value, say 1e-5. Make sure the case diverges quickly. You should see your monitor points oscillating and diverging. Increase the coefficient an order of magnitude at a time until the oscillations stabilize, then fine tune it to get a critically damped response. Ideally your monitor points should overshoot just slightly, then converge WITHIN each timestep.
Keep us posted.

[Zymon]

Wow, thanks a lot for all your advices. I will try it today and tell you later on about the progress. Another short question: What should I set for the 'Variables (for sources and relevant sink options only) -> Velocity' in the boundary source options? Should I use velocity values of zero or set it to 'Determine from mass fluxes'?

Here are some comments to your post:
- The membrane is fixed at its end and is also fixed to a ring, which also has a fixed support. Furthermore, it has contact with the pressure plate, which pushes the membrane into the chamber. Is there anything else I should consider in the mechanical part of the simulation?
- Yes, right. The structure is not pre-stressed. The actual initial conditions in the chamber are much higher than 1 mm Hg. I reduced the value for testing to prevent a sudden pressure load. The standalone FEM simulation can handle a sudden pressure load of about 20 mm Hg without problems. So this should not be a problem. I also reduced it to 0 mm Hg for testing, no differences. When I get rid of the instabilities, I will use a ramping of the outlet boundary condition. Or is it better to set the the reference pressure to the high value and the outlet pressure to 0?
- Yes, 40 coupling iterations are way too much. You are right. I set it to this high value when I tested a really low under relaxation factor. But this was not a good option, as you know. However, it is advisable to use a minimal number of coupling iterations, right?

I'll keep you updated!

[stumpy]

- Conceptually the mechanical side sounds OK. The constraints on the edge of the membrane sound consistent between the fluid and structural sides.
- If standalone FEM can handle 20 mm Hg (transient, same timestep size, without substepping) then FSI should be OK too, BUT it would potentially need a lot of damping (large coefficient value) to remain stable. After the initial bump that coefficient would probably be too large, so you'd have to stop, reduce the coefficient, then restart. Ramping up the outlet BC from 0 to 20 mm Hg is one option (try to use a large timestep during the ramping so you're not resolving too much of the transient), the other is to do a steady state 2-way FSI to get to the pre-stressed state then re-start this as a transient. That restart process is not obvious.
- A minimal number of Coupling Iterations is advisable. With an URF=1 then at least 2 loops is still a good idea. In theory the convergence criteria should take care of not stopping the coupling iterations until the loads are converged, but it's always good to check that with monitor points.

[Zymon]

Hi stumpy,

just wanted to inform you that it is working now thanks to all of your advices. I also changed the mesh and contact parameters of the FEM setting a bit. Now the structural part is also performing better. Loading the pressure from 1 mm Hg to over 100 mm Hg is also working good. Thanks a lot for your help!

Is there any literature about the mass flux pressure coefficient? I found the papers by Degroote about the artificial compressibility. Is this a comparable approach?

Best wishes!

[stumpy]

I don't know of any papers. Keep in mind that the mass flux pressure coefficient is just the hook that is used in CFX to get access the the [A] matrix coefficients on the continuity equation (in [A][x] = [b]). So any literature would probably not talk about this name. A good starting point in the literature is:
http://www.sciencedirect.com/science...45794908001466
Artificial compressibility sounds similar, but it does depend on how it is done. When you modify a coefficient in [A] you are not changing the converged solution, only the convergence path/stability.

[Zymon]

Can you please explain this in more detail? How does the mass flux pressure coefficent change the convergence path in one time step?

As far as I understand, the instabilities occur because the flow can't respond to the change of the structure due to the incompressibility of the fluid. So the forces at the interface increase tremendously. Changing the MFPC changes the coefficient directly in the [A] matrix of the continuity equation at the interface. And MFPC is the derivative of mass flux with respect to pressure. With higher MFPC the mass transfer from the interface to the solid is damped for each fluid iteration and therefore the pressure response from the fluid part is damped. Is this correct?