[Hiroshiman]

Hi,
I'm running OpenFOAM on cases with backflow occuring at outlets at some times ; it is a physiological flow and therefore cannot be avoided and I cannot change the geometry adding outlet extrusions.
On outlets the pressure is imposed (thanks Groovy) and the zeroGradient BC is used for velocity. When backflow occurs, OpenFOAM offers BCs like pressureNormalInletOutletVelocity but it's neither enough nor desirable to have a normal velocity entering the domain.

One idea and what I'm trying to do :
if flow leaves the domains
then zero traction
else substract a traction term depending on the normal component of the in-going velocity, local to the considered face.
The adding term is an outward traction pushing the flow in the outward normal direction.
In grooveBC, I tried something like :

type groovyBC;
variables"n=normal()/mag(normal());Unn=U&normal()*n;K=10;";
gradientExpression      "((phi>0) ? vector(0,0,0) : Unn)";
fractionExpression "0";
value uniform (0 0 0);

I tried for different values of K and always end up diverging (most of the time later than if I were using only zeroGradient). Does any of you have an idea on how to improve this ?

Best regard,
Florian

[Hiroshiman]

Any swak users here have an idea on how to impose a custom traction BC ?

[gschaider]

Originally Posted by Hiroshiman

Hi,
I'm running OpenFOAM on cases with backflow occuring at outlets at some times ; it is a physiological flow and therefore cannot be avoided and I cannot change the geometry adding outlet extrusions.
On outlets the pressure is imposed (thanks Groovy) and the zeroGradient BC is used for velocity. When backflow occurs, OpenFOAM offers BCs like pressureNormalInletOutletVelocity but it's neither enough nor desirable to have a normal velocity entering the domain.

One idea and what I'm trying to do :
if flow leaves the domains
then zero traction
else substract a traction term depending on the normal component of the in-going velocity, local to the considered face.
The adding term is an outward traction pushing the flow in the outward normal direction.
In grooveBC, I tried something like :

type groovyBC;
variables"n=normal()/mag(normal());Unn=U&normal()*n;K=10;";
gradientExpression      "((phi>0) ? vector(0,0,0) : Unn)";
fractionExpression "0";
value uniform (0 0 0);

I tried for different values of K and always end up diverging (most of the time later than if I were using only zeroGradient). Does any of you have an idea on how to improve this ?

I feel a bit stupid when looking at the equation: I don't see where K is used so it doesn't surprise me that it has no effect.

I'm not sure whether an expression like this would help (gradients are really tricky and you give a zero-gradient in non-normal-direction which basically means "do what you want"). What might be better is giving some kind of dirichlet-condition where flow comes in (the way inletOutlet does)

BTW: normal() is already normalized

[Hiroshiman]

Hi Gschaider,
OpenFOAM offers mixed boundary conditions like pressureNormalInletOutletVelocity but I wanted to add some kid of gradient penalisation on the backflow; Dirichlet conditions seem a bit too invasive for that. The term is to be backflow velocity dependant.
I'll switch to a Dirichlet condition in case of backflow, I'm doing some tests constraining the the profile to be parabolic when entering the domain with groovyBC.
And yes I forgot to multiply the corrective term by "K", my bad.
Thank you for your answers and groovyBC that makes it easy to create unnatural conditions.
Regards.

[gschaider]

Originally Posted by Hiroshiman

Hi Gschaider,
OpenFOAM offers mixed boundary conditions like pressureNormalInletOutletVelocity but I wanted to add some kid of gradient penalisation on the backflow; Dirichlet conditions seem a bit too invasive for that. The term is to be backflow velocity dependant.
I'll switch to a Dirichlet condition in case of backflow, I'm doing some tests constraining the the profile to be parabolic when entering the domain with groovyBC.
And yes I forgot to multiply the corrective term by "K", my bad.
Thank you for your answers and groovyBC that makes it easy to create unnatural conditions.
Regards.

OK. Did the K improve things?

Anyway. It is quite hard to do stabilisation stuff "over the net" (usually the way how it diverges gives an idea how to stabilize it) so you're on your own here. Just one idea: Have a look at storedVariables. What you could do is store the penalisation (Dirichlet or Neuman) from the last time-step and when you calculate a new one use this to "underrelax" it (basically mix old and new one). That might improve things as the correction won't "overreact" (which might lead to oscillations that might diverge)

[Hiroshiman]

Hi,
I ran some tests with various values of K and now the solution is pretty stable, even on a complex case (artery with one inlet and 8 outlets) with an important backflow.
I think I should make K a function of the normal to the wall velocity of the fluid incoming to give a more proportional correction.
I haven't noticed any oscillation yet but I keep with me the idea of using storedVariable for underrelaxing the solution, thank you for the idea.

Best regards,
Florian

[90nash]

Greetings Hiroshiman,

I am trying to simulate flow in a simple lung bifurcation (one inlet and two outlet with outlets of equal radii) with a parabolic velocity type profile at the inlet (fluctuating with a Womersley number of 4.3). I am using pressureNormalInletOutletVelocity BC at both the outlets. However the solution seems to diverge everytime (no matter what i use in the fvSchemes dictionary or even on employing relaxation factors). The solution diverges after the flow reverses at the outlets i.e. after time period of T/2.

I have been stuck on this problem for over a month and was playing around with different schemes. However your post has given me hope. Can you please explain the boundary condition that you have used and why using pressureNormalInletOutletVelocity does not work. I dont quite understand it. It seems using this BC seems to get you to a converged solution. I would really appreciate your help on the matter.

Many Thanks in Advance.
Best,
Nadish

[Hiroshiman]

Hi 90nash,
sorry for the late reply.
First, pressureNormalInletOutletVelocity bc works quite well as long as the flow reversal happens where the flow isn't too perturbed. Typically adding extrusions at your outlets will make the simulation stable most of the time.
If for some reason you do not want to use extrusion, there are few solutions. They are the least thing to do but can help. Read this article:
A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations DOI: 10.1007/s00466-011-0599-0
The proposed condition is derived from an article from Bazilevs and conceived from and for finite-element but I tried to adapt the idea here in OpenFOAM.
The concept is to use groovyBC to apply a zeroGradient bc whenever the flow exit the domain and a Neumann bc when it enters the domain. The idea it to give convective "information" to the entering flow, the gradient value being a function of the incoming velocity normal to the cell. It is local to every boundary faces.
This is not complex but a personal interpretation of the article, use at your own risk.
Florian

[90nash]

Hello Hiroshiman,

Thank you for your reply. The article that you provided is highly illuminating! Apparently divergence due to backflow is quite a common problem. I am trying to set up a case by employing stabilization at the outlet as you did. I will share the results soon.

Meanwhile i had created a mesh with extended outlets.Using pressureNormalInletOutletVelocity at the extended outlets did solve my problem. The only problem that i am facing now is at the inlet. After the flow reversal (exhalation phase, time period > 0.5) near the inlet i am getting bizarre velocity vectors. Similar to what i was getting at outlets before extending the geometry. Please see attached contour plot in which you can notice a "kink" near the inlet patch. I have also attached an image of how vectors look like. This region is far from my area of study (near bifurcation) so i may ignore this as long as it does not affect the result.

Can you guess what might be causing this "kink" at inlet and how to get rid of this?
Will extending the inlet further help remove this "kink"?

I am attaching the inlet velocity BC's i used:

    INLET
    {
        type            groovyBC;
        variables       "Uavg=0.343440945;yp=pos().y;zp=pos().z;R=max(zp);a=2*Uavg/(R*R);vel=((2*Uavg)-(a*((yp*yp)+(zp*zp))))*sin(3.14159265359*time());";// Time period=3 sec
        valueExpression "vector(vel, 0, 0)";
    }

    WALL
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    OUTLET1                                            
    {
        type            pressureNormalInletOutletVelocity;
        phi             phi;
        rho             rho;
        value           uniform (0 0 0);
    }

    OUTLET2
    {
        type            pressureNormalInletOutletVelocity;
        phi             phi;
        rho             rho;
        value           uniform (0 0 0);
    }

[Hiroshiman]

Your inlet velocity is set on the direction x ? Shouldn't it be on z (or face normal) ?

[90nash]

The inlet patch lies in the y-z plane and velocity is in the +x. I guess the coordinate at the bottom of image is confusing because of orientation. But if you notice carefully the patch is indeed in y-z.