Variational Form Navier-Stokes for Windkessel afterload
5 months ago by
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.
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.
link documentation: http://hplgit.github.io/fenics-mixed/doc/pub/fenics-mixed-4paper.pdf
link to scrips: https://github.com/hplgit/fenics-mixed/tree/master/doc/src/src-fenics-mixed/ns
Normal Variational form of Tentative velocity step:
F1 = \ rho*dot((u - u_n) / k, v)*dx \ + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \ + rho*inner(grad(U), grad(v))*dx \ - inner(p*Identity(len(u)), grad(v))*dx \ + dot(p_n*n, v)*ds \ - dot(mu*nabla_grad(U)*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(grad(u_n)*u_n, v)*dx \ + mu*inner(grad(U), grad(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
By dropping this term, one assumes the Neumann BC stating: Pressure in normal direction is zero at outlets and wall. Is this valid?
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)
- dot(mu*nabla_grad(U)*n, v)*ds
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,
Community: FEniCS Project
8 weeks ago by
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!!
Please login to add an answer/comment or follow this question.