Variational Form Navier-Stokes for Windkessel afterload

360
views
0
10 months ago by
Dear Community,

I am trying to solve a CFD problem for a geometry with one intlet and multiple outlets (geometry of the aortic arch). As boundary conditions, I want to prescribe the velocity at the inlet and variable pressure at every outlet. The pressure is calculated for every outlet with a two-element windkessel model as a function of the previous pressure and flow.

Equation:
p(n+1) = p(n) + dt/2*(q(n) - p(n)/R )/C
with p(n+1) the new pressure, p(n) current pressure, q(n) current flow, dt = time step size, R = resistance (constant), C = compliance (constant)

After many attempts, I found some working code that did exactly what I want. (links below)
But I noticed that some terms were missing in the variational form of the NS equation in the tentative velocity step.

Normal Variational form of Tentative velocity step:
F1 = \
rho*dot((u - u_n) / k, v)*dx               \
+ dot(p_n*n, v)*ds                         \
- rho*dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)​

Variational form of Tentative velocity step in working scrips:
F1 = rho*inner((u - u_n)/k, v)*dx \
- rho* inner(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
​

I would like to know what the effect of the missing terms is and what effect it has on the solution.

The following terms are missing:
- inner(p*Identity(len(u)), grad(v))*dx​
effect?

dot(p_n*n, v)*ds​
By dropping this term, one assumes the Neumann BC stating: Pressure in normal direction is zero at outlets and wall. Is this valid?

- dot(mu*nabla_grad(U)*n, v)*ds​
By dropping this term, one assumes the Neumann BC stating: gradient of velocity is zero in normal direction. (i.e. fully developped flow at outlets)
So neglecting this term is legid in my case.

But is it correct to omit the other two terms?

I hope someone can be of help,

Regards,
Martijn
Community: FEniCS Project
The term
- inner(p*Identity(len(u)), grad(v))​​

can be interpreted as the virtual power of pressure. Omitting it effectively neglects a hydrostatic pressure term in the stress tensor.

But your first code snippet should not work anyway since in the line

 - inner(p*Identity(len(u)), grad(v))

corresponding to the pressure you are missing the integration:

 - inner(p*Identity(len(u)), grad(v))*dx

written 10 months ago by klunkean

Thank you for the response. I indeed forgot the intergration term here, but it is just a typing error when writing this question. It is present in the script. (corrected it in the question above)

It seems to me that it is not legit to omit this term for 3D channel flow problem, since the driving force of the pressure is neglected. Is that right?

written 10 months ago by Martijn
If you are referring to the linked file NSSolver.py, then they are not forgetting the pressure but rather considering the pressure to be an independent variable which is updated. The idea is outlined here.

Omitting the term altogether and not having a different form from which the pressure is updated, however, should give unsatisfying results.
written 10 months ago by klunkean

0
7 months ago by
The pressure term in the Navier-Stokes is the contribution to the velocity field due to the change of pressure in some direction. This can also be seen from the fact that this term is expressed as the divergence of the pressure.  In particular, in an x-direction oriented flow the derivative of the pressure with respect to x is not zero. Actually it is negative and as a result the pressure term is not zero.

However, even if you have not taken pressure in consideration in the first step of the tentative velocity, you add it on the third step, which is the same with a small price of accuracy of course. And this is called the Chorin's projection method.

Now, what happens here, is that you have a code using Chorin's method, which solves the tentative velocity step without taking pressure at the last time-step, but deletes it from the equation. That is why your code is working. The scheme that you declare as "Normal Variational form of Tentative velocity step"  is calles IPCS and takes into account the pressure at the last time-step. This scheme is much more efficient than Chorin's.

Now, if you look more closely, you will see in the NSSolver.py that the second step that calculates the pressure is different than the second step of the IPCS scheme that you use to compare. And there lies the missing pressure from the first step. This is where the price is paid for omitting the pressure from our equation. It should be somewhere right?

So you have two different numerical schemes (both working correctly), Chorin's and IPCS and the difference between them is the way they treat pressure. If you try harder you will see that also the IPCS scheme works just fine (and maybe better) for you!!