Boundary term appears to be missing in fenics Navier Stokes example (?)

[mekong]

Dear all,
In this example on a numerical solution of the Navier-Stokes equations with Fenics, a boundary term appears to be missing (?)

Consider the second step in the section ‘Variational formulation’, and the second step in the incremental pressure correction scheme (IPCS). One obtains the equation

-\nabla \cdot u^{\star} / \Delta t + \nabla^2 p^{n+1} - \nabla^2 p^n = 0

then one multiplies by a test function q, integrates, and obtains

-\langle \nabla \cdot u^{\star}, q \rangle_\Omega / \Delta t + \langle \nabla^2 p^{n+1} - \nabla^2 p^n , q \rangle_\Omega = 0, where the subscript \Omega denotes integration within the manifold \Omega.
The second term is then integrated by parts
\langle \nabla^2 p^{n+1} -\nabla^2 p^n , q \rangle_\Omega = - \langle \nabla (p^{n+1} - p^n) , \nabla q \rangle_\Omega + \langle q [\nabla (p^{n+1} - p^n) ] \cdot \hat n \rangle_{\partial \Omega},
where \partial \Omega is the boundary of \Omega and \hat n its unit normal vector pointing outside \Omega.

Thus the variational problem to solve should be

-\langle \nabla \cdot u^{\star}, q \rangle_\Omega / \Delta t - \langle \nabla (p^{n+1} - \nabla) p^n , \nabla q \rangle_\Omega + \langle q\nabla[ (p^{n+1} - p^n) ] \cdot \hat n \rangle_{\partial \Omega} = 0

This is the same equation as (3.34) here, modulo the boundary term, which is missing in there. Why? I don’t see any reason why that term would vanish.

Thank you
Best,

[Mingming_Zhang]

Hi,

So the last term in your last equation is basically grad(p).n = 0. Take right wall of a box as an example, its unit normal vector is (1,0), so the gradient of pressure in x direction must be 0 if grad(p).n = 0 holds. That being said, p is a constant in x direction on the right wall, which means the total force acting upon right wall is balanced (inside and outside of the right wall). So the box is rigid, not deformable, as I interpreted it. Hope it makes sense to you.

[dokken]

Ive written up notes summarizing bcs on splitting schemes at
https://computationalphysiology.github.io/oasisx/splitting_schemes.html

[mekong]

Hello,
I am sorry, but your reply is totally unclear to me, see below

@Mingming_Zhang

So the last term in your last equation is basically grad(p).n = 0.

It is not: the test function appears under the gradient operator too.

@Mingming_Zhang

p is a constant in x direction on the right wall

The right wall is a line joining the two points (L,0) and (L,h), and it has x=const so it does not make sense to state that a field is ‘constant in the x directon on the right wall’.

[Mingming_Zhang]

I’m sorry, is it \nabla (pq) \cdot \vec{n} or q \nabla p \cdot \vec{n}? I guess it should be the latter.

If it is \nabla p \cdot \vec{n} = 0 with \vec{n} = (1, 0), we have \frac{\partial p}{\partial x} = 0 at x=1 (right wall). This is telling you p at both sides of the right wall are identical so that p has no gradient along x direction at x=1, though p might be a function of y. For example,

10.001  | 10.001
9.021   | 9.021
-0.203  | -0.203
2.189   | 2.189

where vertical line | is the right wall.
Forget about it if this does not make sense to you.

[dokken]

@mekong

The second term is then integrated by parts
\langle \nabla^2 p^{n+1} -\nabla^2 p^n , q \rangle_\Omega = - \langle \nabla (p^{n+1} - \nabla) p^n , \nabla q \rangle_\Omega + \langle \nabla[ (p^{n+1} - p^n) q ] \cdot \hat n \rangle_{\partial \Omega}⟨∇2pn+1−∇2pn,q⟩Ω=−⟨∇(pn+1−∇)pn,∇q⟩Ω+⟨∇[(pn+1−pn)q]⋅^n⟩∂Ω\langle \nabla^2 p^{n+1} -\nabla^2 p^n , q \rangle_\Omega = - \langle \nabla (p^{n+1} - \nabla) p^n , \nabla q \rangle_\Omega + \langle \nabla[ (p^{n+1} - p^n) q ] \cdot \hat n \rangle_{\partial \Omega},

Integration by parts here is wrong.

[mekong]

You are right, I corrected my original post.

My original question still stands: q in the boundary term \langle q [\nabla (p^{n+1} - p^n )] \cdot \hat{n}\rangle_{\partial \Omega} vanishes on the parts of \partial \Omega where p is known, i.e., on the right wall. However, in general it does not vanish elsewhere, so I don’t see why this term is missing as a whole in the Fenics example.

[dokken]

Please read through the notes i referenced earlier

@dokken

Ive written up notes summarizing bcs on splitting schemes at
https://computationalphysiology.github.io/oasisx/splitting_schemes.html

[Mingming_Zhang]

I guess this is not a Dirichlet BC where the field is prescribed (or is given numbers directly) so that test function vanishes. In addition at least p^{n+1} is not known.

[mekong]

I did, but I don’t see the answer in there.

[dokken]

See specifically:
https://computationalphysiology.github.io/oasisx/splitting_schemes.html#essential-boundary-conditions
and the following section.
Note that \phi ln those notes represent the difference between the old and new pressure

[mekong]

@dokken

See specifically:
https://computationalphysiology.github.io/oasisx/splitting_schemes.html#essential-boundary-conditions
and the following section.

I could not find the answer in there. It says “However, if \partial \Omega_N \neq 0, then we need to consider the open boundary conditions for the pressure correction schemes”. In the channel flow example of my original post \partial \Omega \neq 0, so we are in this case.

However, still I don’t see why the boundary terms are absent in the code for the channel flow of my original post. I wonder whether this is because in the pressure-correction scheme, one assumes that p^{n+1} takes a known value on \partial \Omega, and thus q= 0 on \partial \Omega (?)

[dokken]

The text states that at all Dirichlet conditions for velocity, one tend to use dphi/dn=0 on this boundary.

For the remaining boundary the one where we have outflow, you set phi=0. Thus you have no boundary term.

[Mingming_Zhang]

So what is the physical meaning of \nabla (p^{n+1} - p^n) \cdot \vec{n} = 0? pressure gradient normal to the wall stays the same? Could you show us an example? Thanks!

[mekong]

Thank you for your reply. Of course, if \partial \phi / \partial n= 0 on the boundaries where one has Dirichlet conditions for the velocity, then that boundary term vanishes.

@dokken

The text states that at all Dirichlet conditions for velocity, one tend to use dphi/dn=0 on this boundary

That is not stated in the screenshot that you posted.

For the sake of other users, and in an effort to improve this beautiful Fenics project, I would like to point out some clarity issues on this point:

• The only place in this webpage where the condition \partial \phi / \partial n= 0 on \partial \Omega is stated is in Section “Essential boundary conditions,” paragraph “Assume that \partial \Omega_N = 0 […] Thus we use that \partial \phi / \partial n = 0 on \partial \Omega.” I find this very unclear, because, as it is, the text implies that the condition \partial \phi / \partial n = 0 on \partial \Omega holds only within the conditions of that paragraph, i.e., in the case where \partial \Omega_N = 0 . Also, “Thus we use that \partial \phi / \partial n = 0 on \partial \Omega” does not specify that this condition holds only for boundaries where Dirichlet boundary conditions for the velocity are imposed.
• Something about this boundary term should be said in the fenics example program or webpage on Navier Stokes equations, it is very confusing to have to go through all this to find out eventually that that tern vanishes because an unmentioned boundary condition was hypothesized.

I hope that this will help!

[dokken]

@mekong

That is not stated in the screenshot that you posted.

That is stated in the text i first referenced to, that references the second screenshot

As you refer to.

If you continue looking through this text, you will get the full formulation.
I.e. The equations in:
https://computationalphysiology.github.io/oasisx/splitting_schemes.html#navier-stokes-equation

Im open for suggestions on how to improve the text though.

@mekong

Something about this boundary term should be said in the fenics example program or webpage on Navier Stokes equations, it is very confusing to have to go through all this to find out eventually that that tern vanishes because an unmentioned boundary condition was hypothesized.

Please note that you are looking at an 8 year old tutorial, while the webpage i am referencing is maintained (by me).

I Also have an updated version of the splitting scheme at:
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html
which writes out the assumptions on the boundary condition explicitly.