[newOpenfoamUser]

Hi,

I am running laminar simpleFoam simulations (2 orthogonal correctors) on geometry of blood vessels. The vessels typically have 1 or 2 inlets and multiple outlets. I am simulating the cases with different inlet velocities. When the velocity is high, I observed some cases start to have reverse flow at the outlet and the simulation slowly diverged. Whereas it almost doesn't happen when velocity are lower. There are 3 reasons of divergence that I am suspecting now.

The first is that the velocity is too high and there might be some turbulence and the laminar static simulation cannot capture the effect. However, I am not sure whether this will make it diverge. I calculated the Reynold numbers for a few cases with high inlet velocity, based on the formula Re=(Diameter of inlet)*(maximum velocity at inlet)/(kinematic viscosity). In 2 cases, I found the Reynold number to be 2000 and 2200, seems like close to onset of turbulence.

The second reason I suspect is the boundary condition I set is causing the reverse flow and divergence. The boundary condition that I use is parabolic velocity at inlets, zeroGradient boundary condition for pressure at inlets. At the outlets, I use zeroGradient for velocity. For the pressure outlet condition, I use a special case of windkessel model (https://en.wikipedia.org/wiki/Windkessel_effect). Basically, I set C=0 in the 2 element windkessel model. Therefore, the pressure at each outlets are given by P=QR, where R is a prespecified resistance value, Q is the volumetric flow rate at each outlets. I am using groovyBC to implement the pressure boundary conditions at outlet. It looks like this

outlet1
    {
	type            groovyBC;
	variables	"res=10000;flow=sum(U&Sf());";
	valueExpression	"flow*res";
    }

I do not know whether this implementation is correct and hope someone can point it out if this is not correct. I am not very familiar with the details of FVM. After browsing the forum, I read that there are implicit and explicit way of implementing boundary conditions and implicit way is more stable. So, I am wondering whether this is the case where it should be implement implicitly and I have done it explicitly?

The third reason might be that I do not extend the domain enough. While this is true for some case, for some cases it does seems like it is extended enough. See the figure attached below. Vessel1 have 2 outlets that are close to a bifurcation whereas the outlets of vessel2 are far from the bifurcation.

I think the mesh should not be the problem here. For both meshes, the number of cells are about 300k. I attached the checkMesh output below. The output seems to suggest they are ok.

Hope someone can help me with this. Thanks!

[newOpenfoamUser]

Something that I forgot to mention is that I have tried using inletOutlet for velocity outlet boundary condition. This seems like a common thing to do in the case of reversed flow at outlets. It kind of stabilised the simulation and prevent it from diverge. However, I have a feeling that this does not really solve the problem and there is still some fundamental issues. In particular, if I move a little bit inwards from the outlet, there is still some area with the velocity pointing inwards. Also, the residuals do not seem to decrease to small value. It might stabilised at some relatively large value like 0.01 (for p) and 0.001 (for U).

[piu58]

Zero gradient is be the right b.c. for outlets *in principle*. But it may be, that you set this condition to close to the turbulencey your model calculates. I recommend enlarging the outflow area with a tube of constant diameter. In this tube, the flow has time to settle.

In effect, you cannot get a zero gradient condition this near to disturbations of the flow. This is nonphysical. It may be apllicable with low Re numbers, but starts making problems at turbulent flow.