[opencurrent-users] Boundary Conditions for Navier-Stokes Equation

[Ben Mueller]

I am having trouble in changing the boundary conditions of the fluid
flow such that the flow is not limited to a closed box. I am
simulating a smoke plume which should rise and exit the simulation
domain. I am basing my code on the cubicrayleigh.cpp example code.

Is there a general explanation of all the different boundary condition
types and how to apply them correctly?

[Jonathan Cohen]

Implementing open boundary conditions shouldn't be too bad, although it is somewhat different from what is currently there.

the basic idea can be found in this paper:

http://jcohen.name/papers/Shah_Extended_2004.pdf

section 3.2

as far as I know, this isn't well described in the graphics literature.  Non-graphics people do stuff like this all the time, but they usually use something more sophisticated.

The basic idea is to enforce du / dn = 0 at all boundaries.  This can be handled in the ghost cells by setting the ghost cell value based on the derivative of the interior cell value during the apply_3d_boundary_conditions function, e.g.:

u[0] - u[-1]  = u[1] - u[0]

=> u[-1] = u[1] - 2 * u[0]

you then need to enforce this boundary condition during the multigrid iterations.

there are two approaches to this.  The easy approach is to pretend that the boundaries are simply wall boundaries, and then reset the ghost cells based on the du/dn=0 condition at the end.  This will result in some divergence along the border cells, but basically works.

The 'correct' approach is to modify the relaxation step slightly. 
The scheme is slightly obfuscated, but what is currently done is described here:
http://code.google.com/p/opencurrent/wiki/EnforcingBoundaryConditionsInMultigrid

Basically, you can come up with a different set of weights in the relax step based on dU/dn = 0, and modify relaxation for the border cells appropriately.

Let me know if this makes sense, or if you need more clarification.

[Jonathan Cohen]

sorry, i was being dense before.  here's what you should actually do:

The correct boundary conditions are simply: u[0] = u[1]

that is, you set the flux on the external face to be equal to the flux one grid cell in.  And then you probably want no-slip for the tangential components.  And since the ghost-cells need to be filled in, set u[-1] = u[0].

Then, when you solve for pressure and subtract grad p, it will overwrite u[1], so that u[0] no longer equals u[1].  So you have a choice - either reinforce the 'open' boundaries which will result in some non-zero divergence for border cells, or leave things alone which results in zero divergence but the du/dn=0 condition not being satisfied.

In my experience, it is usually better to have div u = 0 than accurate boundary conditions.

And then the 'correct' but more involved solution is to modify the Poisson solver to solve for P while also enforcing that d^2P/(dn)^2 = 0, which would satisfy both conditions when grad p is subtracted.  but that is probably more trouble than it's worth for graphics.

[Ben Mueller]

thank you very much for your help.

unfortunately I don't get the correct results.

I set the bc's to BC_FORCED_INFLOW_VARIABLE_SLIP with value=0.0 and
aux_value=0.0 and then adjust the bc's for the top of my simulation
domain just before the call to the projectionSolver.
I would set v[i,ny-1,k] = v[i,ny-2,k] and v[i,ny,k] = v[i,ny-1,k]
(ghost cells) for all i,k.
for the no-slip of the tangential components I would set u[i, ny-1,k]
= 0 and w[i,ny-1,k] = 0, so there should be no flow in the tangential
direction.

the flow is moving upwards and then 'hits' the ceiling and moves
tangentially to the ceiling.

what am I missing?

[Jonathan Cohen]

in the y direction, i think you want

v[ny] = v[ny-1]
v[ny+1] = v[ny-1]

since there are ny+1 cells in that direction (staggered grid).