Boundary conditions in Navier Stokes demo

Boundary conditions in Navier-Stokes 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://onlinelibrary.wiley.com/doi/10.1002/fld.1650071008/abstract>.

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.

On 8 April 2013 00:55, Peter Franek <email address hidden>wrote:

> New question #226141 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/226141
>
> In the incompressible Navier-Stokes 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
> Navier-Stokes 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.
>

On Sun, 07 Apr 2013 23:46:21 -0000
Nico Schlömer <email address hidden> wrote:
> Question #226141 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/226141
>
> Status: Open => Answered
>
> Nico Schlömer proposed the following answer:
> Boundary conditions in Navier-Stokes 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 non-zero 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://onlinelibrary.wiley.com/doi/10.1002/fld.1650071008/abstract>.
>
> 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.

On 9 April 2013 18:41, Jan Blechta <email address hidden>wrote:

> Question #226141 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/226141
>
> 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://answers.launchpad.net/fenics/+question/226141
> >
> > Status: Open => Answered
> >
> > Nico Schlömer proposed the following answer:
> > Boundary conditions in Navier-Stokes 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 non-zero 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://onlinelibrary.wiley.com/doi/10.1002/fld.1650071008/abstract>.
> >
> > 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), 1125-1146, 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/projection.

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.
>

Thank you guys for the comments. Which method would you recommend me for incompressible flow with an everywhere non-slip 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 divergence-free 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

On 15 April 2013 01:11, Peter Franek
<email address hidden>wrote:

> Question #226141 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/226141
>
> 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 non-slip 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
fenics-book.
http://link.springer.com/chapter/10.1007/978-3-642-23099-8_21
Here, we tested it for the lid-driven cavity problem. (Note that you
can use the normalize function. )

>
> Another question: can I not define a functionspace of divergence-free
> 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 divergence-free 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.
>

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=UnitSquareMesh(32,32)
f = Expression(("1", "0"))

# Define function space and basis functions
V = VectorFunctionSpace(mesh, "CG", 2)
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*inner(grad(u), grad(v))*dx - inner(f, v)*dx
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(grad(p_new), v)*dx

# 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/velocity_nv.pvd")
pfile = File("output/pressure_nv.pvd")

# Time-stepping
t = dt
while t < T + DOLFIN_EPS:

    # Update pressure boundary condition
    #p_in.t = t

    # Compute tentative velocity step
    begin("Computing tentative velocity")
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u1.vector(), b1, "gmres", "default")
    end()

    # Pressure correction
    begin("Computing pressure correction")
    b2 = assemble(L2)
    solve(A2, p1.vector(), b2, "gmres", prec)
    (p_new,c)=p1.split()
    end()

    # Velocity correction
    begin("Computing velocity correction")
    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()
-------------------------------------

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Please, could you give me some hint or reference on how to efficiently solve NS with an everywhere non-slip 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.

On Sun, 05 May 2013 12:11:32 -0000
Peter Franek <email address hidden> wrote:
> Question #226141 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/226141
>
> Peter Franek posted a new comment:
> Please, could you give me some hint or reference on how to efficiently
> solve NS with an everywhere non-slip 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 neumann-poisson demo.

2. fix pressure in one vertex:
   use DirichletBC(pressure_space, 0.0, pin_point, "pointwise") where
   pin_point is SubDomain consisting of one mesh vertex. This has a
   good meaning for CG1 pressure_space...

Jan