Boundary conditions in Navier Stokes demo
In the incompressible NavierStokes demo, you use the Chorin projection method: first you calculate the tentative velocity, then the pressure update. The pressure update is a Poisson equation. This equation, in the variational formulation, yields a unique pressure p_n such that p_n satisfies the Dirichlet boundary condition on [bcp], as well as a von Neumann condition on zero normal derivative in the rest of the boundary of the domain [bcu]. Is it right?
My question is the following: can I prescribe a Dirichlet boundary condition on one part of the boundary [bcu] for the velocity, a Dirichlet boundary condition for the pressure on another part of the boundary [bcp], and FURTHER require a Neuman boundary condition for the pressure on the first part of the boundary [bcu]? Isn't this too much, does such NavierStokes problem have a solution? Peter
Question information
 Language:
 English Edit question
 Status:
 Answered
 Assignee:
 No assignee Edit question
 Last query:
 Last reply:
Revision history for this message

#1 
Boundary conditions in NavierStokes are a tricky business. There is general agreement that the velocity should be zero on the boundary (whenever there's no in or outflow)  in both tangential and normal direction of the boundary that is. What you do with the pressure around the boundaries is actually unclear though.
Check out, for example <http://
When you finally write down the equations, then you have to make a choice, and the way that it is formulated in the demo is:
1. Solve a weak problem for the tentative velocity.
2. Solve a weak problem for the pressure.
3. Solve a weak problem for the update.
All three steps need to be equipped with appropriate boundary conditions.
Important for the velocity are of course the boundary conditions in the last step. In the demo, it is chosen to enforce those with Dirichlet. The price you have to pay for accurate boundary conditions in u is that div(u)=0 is not fulfilled exactly in the vicinity of the boundary.
Revision history for this message

#2 
On 8 April 2013 00:55, Peter Franek <email address hidden>wrote:
> New question #226141 on FEniCS Project:
> https:/
>
> In the incompressible NavierStokes demo, you use the Chorin projection
> method: first you calculate the tentative velocity, then the pressure
> update. The pressure update is a Poisson equation. This equation, in the
> variational formulation, yields a unique pressure p_n such that p_n
> satisfies the Dirichlet boundary condition on [bcp], as well as a von
> Neumann condition on zero normal derivative in the rest of the boundary of
> the domain [bcu]. Is it right?
>
> My question is the following: can I prescribe a Dirichlet boundary
> condition on one part of the boundary [bcu] for the velocity, a Dirichlet
> boundary condition for the pressure on another part of the boundary [bcp],
> and FURTHER require a Neuman boundary condition for the pressure on the
> first part of the boundary [bcu]? Isn't this too much, does such
> NavierStokes problem have a solution? Peter
>
>
Yes, this is exactly what you will do. For splitting schemes you will need
one more bc then the actual equations
permit. Typically, you will have Neuman conditions for pressure where you
have Dirichlet for velocity and vica versa.
But, as you say, this is one condition to many and it may show up as an
unphysical boundary layer.
This boundary layer will typically be O(1) for Chorin, but O(dt) for IPCS
(inc pressure correction scheme),
see also nsbench and the associated study in the FEniCS book.
Kent
> 
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>
Revision history for this message

#3 
On Sun, 07 Apr 2013 23:46:21 0000
Nico Schlömer <email address hidden> wrote:
> Question #226141 on FEniCS Project changed:
> https:/
>
> Status: Open => Answered
>
> Nico Schlömer proposed the following answer:
> Boundary conditions in NavierStokes are a tricky business. There is
> general agreement that the velocity should be zero on the boundary
> (whenever there's no in or outflow)  in both tangential and normal
That's not true. On moving walls, symmetry planes or walls with
Naviers's slip bc etc. there is nonzero tangential velocity.
> direction of the boundary that is. What you do with the pressure
> around the boundaries is actually unclear though.
>
> Check out, for example
> <http://
>
> When you finally write down the equations, then you have to make a
> choice, and the way that it is formulated in the demo is:
>
> 1. Solve a weak problem for the tentative velocity.
> 2. Solve a weak problem for the pressure.
> 3. Solve a weak problem for the update.
>
> All three steps need to be equipped with appropriate boundary
> conditions.
Not necessarily as 3rd step is zeroth order PDE for velocity it does
not need bcs and you may want to keep small error in tangetial velocity
introduced by step 3 with no bcs because of stability and/or
divergence error. Note that typically step 2 is solved solved with
homogenous Neumann so normal velocity is not altered by step 3 with no
bc. See for example
David L. Brown, Ricardo Cortez, Michael L. Minion: Accurate Projection
Methods for the Incompressible Navier–Stokes Equations, Journal of
Computational Physics 168, 464–499 (2001)
> Important for the velocity are of course the boundary
> conditions in the last step. In the demo, it is chosen to enforce
> those with Dirichlet. The price you have to pay for accurate boundary
> conditions in u is that div(u)=0 is not fulfilled exactly in the
> vicinity of the boundary.
Revision history for this message

#4 
On 9 April 2013 18:41, Jan Blechta <email address hidden>wrote:
> Question #226141 on FEniCS Project changed:
> https:/
>
> Jan Blechta proposed the following answer:
> On Sun, 07 Apr 2013 23:46:21 0000
> Nico Schlömer <email address hidden> wrote:
> > Question #226141 on FEniCS Project changed:
> > https:/
> >
> > Status: Open => Answered
> >
> > Nico Schlömer proposed the following answer:
> > Boundary conditions in NavierStokes are a tricky business. There is
> > general agreement that the velocity should be zero on the boundary
> > (whenever there's no in or outflow)  in both tangential and normal
>
> That's not true. On moving walls, symmetry planes or walls with
> Naviers's slip bc etc. there is nonzero tangential velocity.
>
> > direction of the boundary that is. What you do with the pressure
> > around the boundaries is actually unclear though.
> >
> > Check out, for example
> > <http://
> >
> > When you finally write down the equations, then you have to make a
> > choice, and the way that it is formulated in the demo is:
> >
> > 1. Solve a weak problem for the tentative velocity.
> > 2. Solve a weak problem for the pressure.
> > 3. Solve a weak problem for the update.
> >
> > All three steps need to be equipped with appropriate boundary
> > conditions.
>
> Not necessarily as 3rd step is zeroth order PDE for velocity it does
> not need bcs and you may want to keep small error in tangetial velocity
> introduced by step 3 with no bcs because of stability and/or
> divergence error.
Agree
> Note that typically step 2 is solved solved with
> homogenous Neumann so normal velocity is not altered by step 3 with no
> bc. See for example
> David L. Brown, Ricardo Cortez, Michael L. Minion: Accurate Projection
> Methods for the Incompressible Navier–Stokes Equations, Journal of
> Computational Physics 168, 464–499 (2001)
>
Agree, just like to point out that we discussed this to some extent
earlier:
Numerical methods for incompressible viscous flow,
Langtangen, Mardal, Winther,
Advances in Water Resources
25(8), 11251146, 2002
Here, we point out the cause and consequences of the additional and
(slightly) inconsistent
bcs. This is a topic that is often only barely discussed. It seems that
some inconsistency is
inherited by the use of splitting/
Kent
>
> > Important for the velocity are of course the boundary
> > conditions in the last step. In the demo, it is chosen to enforce
> > those with Dirichlet. The price you have to pay for accurate boundary
> > conditions in u is that div(u)=0 is not fulfilled exactly in the
> > vicinity of the boundary.
>
>
> 
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>
Revision history for this message

#5 
Thank you guys for the comments. Which method would you recommend me for incompressible flow with an everywhere nonslip boundary condition? The Chorin method doesn't work (I defined a supplementary normalization condition that the integral of the pressure is zero); even in case of a constant force (the velocity solution should be always zero) it diverges. Do you think that it will go better with IPCS?
Another question: can I not define a functionspace of divergencefree elements? It seems to me that if you have both Trialfunction and Testfunction in such space, the incompressible NS equation simplifies and the pressure terms disappears.
Thanks, Peter
Revision history for this message

#6 
On 15 April 2013 01:11, Peter Franek
<email address hidden>wrote:
> Question #226141 on FEniCS Project changed:
> https:/
>
> Status: Answered => Open
>
> Peter Franek is still having a problem:
> Thank you guys for the comments. Which method would you recommend me for
> incompressible flow with an everywhere nonslip boundary condition? The
> Chorin method doesn't work (I defined a supplementary normalization
> condition that the integral of the pressure is zero); even in case of a
> constant force (the velocity solution should be always zero) it
> diverges. Do you think that it will go better with IPCS?
>
IPCS should be better.
But I would expect Chorin to work also, see e.g. Chap 21 in the
fenicsbook.
http://
Here, we tested it for the liddriven cavity problem. (Note that you
can use the normalize function. )
>
> Another question: can I not define a functionspace of divergencefree
> elements? It seems to me that if you have both Trialfunction and
> Testfunction in such space, the incompressible NS equation simplifies
> and the pressure terms disappears.
>
>
FEniCS does not have any divergencefree elements. The main reason
is that they result in matrices with worse condition number than other
elements. But you can easily set up mixed methods that are weakly
divergence free. See the GRPC method in Chap 21.
Kent
>
> 
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>
Revision history for this message

#7 
Concerning Chorin, I slightly modified the demo and have a noslip bc everywhere, and a constant force f=(1,0). The only reasosonable solution of this problem shoud be trivial (u=0), but it calculates a velocity in the direction of the force and diverges in time. The divergence of the velocity is highly nonzero near the boundary. Could you check please, if I don't have some formal mistake somewhere? I basically changed only the pressure computation step (there is no bc for pressure)
If you don't have time, ignore this please. Anyway, thanks for all your time and for all the answers! Peter

from dolfin import *
k=Constant('0.024')
mesh=UnitSquare
f = Expression(("1", "0"))
# Define function space and basis functions
V = VectorFunctionS
Q = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
P = Q*R
# Define trial and test functions
u = TrialFunction(V)
(p,c) = TrialFunction(P)
v = TestFunction(V)
(q,d) = TestFunction(P)
# Set parameter values
dt = 0.01
T = 3
nu = 0.01
def u_boundary(x, on_boundary):
return on_boundary
nonslip = DirichletBC(V,('0', '0'), u_boundary)
bcu = [nonslip]
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(P)
p_new=Function(Q)
# Define coefficients
k = Constant(dt)
# Tentative velocity step
F1 = (1/k)*inner(u  u0, v)*dx + inner(grad(u0)*u0, v)*dx + \
nu*
a1 = lhs(F1)
L1 = rhs(F1)
# Pressure update
a2 = (inner(grad(p), grad(q)) + c*q + p*d)*dx
L2 = (1/k)*div(u1)*q*dx
# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx  k*inner(
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Use amg preconditioner if available
prec = "default"
# Create files for storing solution
ufile = File("output/
pfile = File("output/
# Timestepping
t = dt
while t < T + DOLFIN_EPS:
# Update pressure boundary condition
#p_in.t = t
# Compute tentative velocity step
begin(
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "default")
end()
# Pressure correction
begin(
b2 = assemble(L2)
solve(A2, p1.vector(), b2, "gmres", prec)
(p_
end()
# Velocity correction
begin(
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "default")
end()
# Plot solution
# plot(p_new, title="Pressure", rescale=True)
plot(u1, title="Velocity", rescale=True)
# plot(div(u1), title="Velocity", rescale=True)
# Save to file
ufile << u1
pfile << p1
# Move to next time step
u0.assign(u1)
t += dt
print "t =", t
interactive()

Revision history for this message

#8 
This question was expired because it remained in the 'Open' state without activity for the last 15 days.
Revision history for this message

#10 
Please, could you give me some hint or reference on how to efficiently solve NS with an everywhere nonslip bc? My problem was that in the Chorin method, the pressure equation has no solution in general, if I use a von Neumann bc for the pressure on the whole boundary. (laplace u = f with von Neumann bc has a solution iff (\int f dx) is equal to (\int (normal derivative of u) dS), or so). I'm trying to modify the equation to include some boundary integrals in the pressure calculation step (so that the above condition is satisfied), but it doesn't behave nicely.
Revision history for this message

#11 
On Sun, 05 May 2013 12:11:32 0000
Peter Franek <email address hidden> wrote:
> Question #226141 on FEniCS Project changed:
> https:/
>
> Peter Franek posted a new comment:
> Please, could you give me some hint or reference on how to efficiently
> solve NS with an everywhere nonslip bc? My problem was that in the
> Chorin method, the pressure equation has no solution in general, if I
> use a von Neumann bc for the pressure on the whole boundary.
> (laplace u = f with von Neumann bc has a solution iff (\int f dx) is
> equal to (\int (normal derivative of u) dS), or so). I'm trying to
> modify the equation to include some boundary integrals in the
> pressure calculation step (so that the above condition is satisfied),
> but it doesn't behave nicely.
>
There are at least 2 approaches:
1. use lagrange multiplier to enforce \int p dx = 0 constraint:
see neumannpoisson demo.
2. fix pressure in one vertex:
use DirichletBC(
pin_point is SubDomain consisting of one mesh vertex. This has a
good meaning for CG1 pressure_space...
Jan
Can you help with this problem?
Provide an answer of your own, or ask Peter Franek for more information if necessary.