Computation of a Simplified Centrifugal Blood Pump

[SRIB]

Hi everybody,
currently I am working in Prof. Behbahani's group on a CFD benchmark study released by the Food and Drug Administration USA that is related to a simplified centrifugal blood pump. The purpose of this project is to determine how computational fluid dynamics can be effectively used to characterize fluid flow and to predict blood damage in medical devices. This pump consists out of two parts, one rotor and one housing. In addition i created a boundary cylinder which surrounds the rotor.
The groups of the stator are inlet, outlet, fixed wall and boundary. The groups of the rotor are wall (the rotor itself, also fixed) and boundary. To get a rotation i set OmegaZ in „Coriolis Source Terms“ to a certain rotation in s-1.
In order to let the flow compute i created the following study:
./code_saturne –s Pump –c Rotor –c Stator
then i set the parameter in the GUI for both rotor and stator. In the subroutine cs_user_coupling.c i set the following for the stator:
/*-------------------------------------------------------------------------
* Example 2: coupling with instance "SATURNE_03".
*
* - coupled faces of groups "coupled_faces"
* - coupled cells (every cell overlapping the distant mesh)
* - all cells available as location support
*-------------------------------------------------------------------------*/
if (true)
cs_sat_coupling_define("Rotor",
"coupled_boundary",
"all[]",
NULL,
"all[]",
verbosity);
}

END_C_DECLS

Fort he rotor i changed the cs_sat_coupling definition to Stator. Checking the listing file the sat/sat coupling seems to work. But after 2 iterations i recieved the following error:

Jacobi [VelocityX]: divergence after 2 iterations:
initial residual: 1.5700e+18; current residual: 3.6045e+23



/home/caelinux/opt/installer/code_saturne-3.0.1/src/alge/cs_sles.c:2583: Fatal error.

Jacobi: error (divergence) solving for VelocityX


Call stack:
1: 0x7f2659e44f1b <reslin_+0x4bb> (libsaturne.so.0)
2: 0x7f2659e641a6 <invers_+0x1f6> (libsaturne.so.0)
3: 0x7f2659c53097 <coditv_+0x1a7b> (libsaturne.so.0)
4: 0x7f2659d329d3 <predvv_+0x56bb> (libsaturne.so.0)
5: 0x7f2659d174ae <navstv_+0x19f6> (libsaturne.so.0)
6: 0x7f2659d53a85 <tridim_+0x45d9> (libsaturne.so.0)
7: 0x7f2659c30014 <caltri_+0x30e0> (libsaturne.so.0)
8: 0x7f2659c06905 <cs_run+0xa35> (libsaturne.so.0)
9: 0x7f2659c05dea <main+0x14a> (libsaturne.so.0)
10: 0x7f265981176d <__libc_start_main+0xed> (libc.so.6)
11: 0x400769 <> (cs_solver)
End of stack

I am reletivly new in code_saturne. So, maybe someone of you has an idea how to fix this problem. Thank you very much for any kind of relevant comments.

[Yvan Fournier]

Hello,

The divergence may be due to one of several possible causes :

- bad mesh quality
- too strong time step
- incorrect setup (you might double-check you have Coriolis and Omega Z only for the rotor, not for the stator)

Did you try running your computation with no rotation ?

A good way of testing the mesh would be to run another computation, joining both meshes in a single domain, with no rotation with a forced inlet velocity, so as have a first time step adjustment (my suggestions, but you may want to check the best practice guidelines on this site for more details: for reasonable CFL, mean should not be above 5 or at the very worst 20, max above 20 might be dangerous; rotation may later require using a smaller time step, depending on rotational velocity)

Finally, you may be interested in out new rotor-stator algorithm, in intermediate version 3.2.1, which we released yesterday:

- rotation is based on mesh joining, updated at each time step.
- this assumes a relatively clean mesh where the rotor and stator meet (ideally, extruded perpendicular to the rotation axis in that area) to avoid joining issues
- no sub-iterations are required, and conservativity between the 2 domains is ensured
- no code coupling required : everything is in one case
- setup assumes using the cs_user_turbomachinery.c user subroutine (simpler in my opinion than cs_coupling.c, but you can check it at http://code-saturne.org/viewvc/saturne/ ... iew=markup to make your own opinion.
- a "frozen-rotor" option also now exists

I am not the expert for this part of the code, but I'll suggest to the specialists to check the forum in case you have more questions.

Regards,

Yvan

[SRIB]

Hello,

thank you very much. Ok, first i will compute it with no rotation in one domain. Does this mean that i have to add both meshes in one case?

The new version is indeed very interesting in terms of my task. I will install it. Finally, what is the frozen-rotor method?

Best regards

[Yvan Fournier]

Hello,

You can either simply use a zero rotation velocity or put everything in one domain, but putting everything in one domain first, and joining the meshes (using non-conforming joining) is better for a first test (and is quite similar to what you would do with the new version anyways).

I believe the frozen rotor method is a method for an approximate "steady solution", so it might not be appropriate, but again, I need to check with the colleagues who implemented this.

Regards,

Yvan

[SRIB]

Hello Yvan,

I`ve more or less three questions:

1.I installled the new version 3.2.1 and wanted to do exactly the same case that i described in the post yesterday but with zero rotational velocity. Although, i used exactly the same subroutine for the coupling it failed. But i checked the subroutine twice and it was exactly the same. I attached the listing file and the subroutine.

2.I let the flow compute with a zero rotation velocity in version 3.0.1 as you mentioned and it stopped after 28 iteations with the following error:

Jacobi: error (divergence) solving for Dissip

Should I now reduce the time step. I used the given time step of 0.1 seconds. I attached the listing file for the rotor. The stator seems to be ok.

3.Yesterday you wrote that I could put everything in one domain and joining the meshes using non-conforming joining. Because of the fact that I am really new in Code Saturne I don´t have an idea how to do this. I checked the SRC folder and found a subroutine called cs_user_mesh. Do I have to use this subroutine and if I have to use it could you perhaps explain how to set it up? In Addition, I don`t really understand the cs_user_turbomachinery.

Thank you very much for your help. This is really very nice.

Regards,

Sebastian

[Yvan Fournier]

Hello,

I didn't check your files yet, but for quick answers:

1)
You won't need cs_user_coupling.c in this version

2) and 3)

Just create one case (either using "code_saturne create" or directly through the GUI), and add the 2 meshes to that case.

For the first test, using the GUI, select both meshes, define a mesh joining, using the selected face groups of both meshes as criteria. Then run the test.

To add rotation, do not define a mesh joining (at least not for the rotor/stator boundary), and define the joining in cs_user_turbomachinery.c instead (also define which cells belong to the rotor before that ; simply adapt the example).

Regards,

Yvan

[SRIB]

Hi Yvan,

i tried to join the meshes. I added both meshes to the case and for the selection criteria i chose the name of my boundary group. It looked like this

Faction: 0.1 Plane: 25 Verbosity: 1 Viisualiation:1 Selection criteria: coupled_boundary

I don`t know if this is right, but i received an error that said that two cells have no connectivity. I read in the user guide about non-conforming mesh joining but at the moment it is hard to understand.

Regards,

Sebastian

[SRIB]

Hi Yvan,

excuse me, please. What I wrote in my third post under point 2 is nonsense, of course. I used a time step (number of iterations) of 50. And after 28 time steps the computation failed. I think this is exactly what you said. Everything above 20 is dangerous.

I did again the same computation with a sat/sat coupling in 3.0.1. I set rotational velocity to 0 and I used 20 time steps. This was successful. Is it possible to try to let the rotor rotate using this case? Which parameter would you then change besides omega z?

I used the following parameter (just for your information):

Rotor:
- unsteady flow
- turbulence model: k-epsilon
- rho=1035 kg/m^3 for blood, viscosity=0.0035Pa.s
- rough wall
-20 iterations
- 4 Processors

next step would be: omega z = 261.799 rad/s (angular pump speed) or 41,6667 s-1 depending on which input code saturne needs.

Stator:

- unsteady flow
- turbulence model: k-epsilon
- rho=1035 kg/m^3 for blood, viscosity=0.0035Pa.s
- Inlet: volumic flow rate: 4.16667E-05 m^3/s, turbulence intensity at inlet: 4%
- rough wall
-20 iterations
- 4 Processors

Regards,

Sebastian

[Jacques Fontaine]

Hello,

Could we have a picture of ur geometry ?

What is the Re number ? Do u really need a turbulence model ?
To start, I think rough walls is not a good idea (you should use smooth walls).

Regards,
JF

[SRIB]

Hello Jaques,

i attached a picture of the geometry.Please let me know if you need other pictures. My reynolds number is 209338.

Regards,

Sebastian