Enforcing a Do-Nothing Boundary

Question

I'm solving a simple 2-D time-dependent navier-stokes problem with an obstruction (circle) located near the inlet. The scheme is crank-nicholson in time, and P2,P1DC (scott vogelius) in space. (Normally, I use a special mesh with these elements, but I can actually reproduce my issue using a delaunay mesh and Taylor Hood elements.) I use Picard iteration for the nonlinear term, which is in rotational form (rot(u)*u).

I have (very) simplified code that reproduces my exact issue below.

My five boundaries are Inflow, Outflow, Top, Bottom, and On the circle. My boundary conditions are no-slip for the circle, top, and bottom boundaries, forced inflow (parabolic), and do-nothing outflow.

It's my (likely incorrect) understanding that I shouldn't have to actually include a term for the do-nothing boundary. However, my results show that my behavior on outflow boundary is strange. If you plot the solution at t=4.0s, you'll see the velocity grows significantly at the center of the outflow. If you force the outflow to have the same condition as the inflow, this problem goes away. (It's easy to add/remove the forced outflow term in the code.)

from dolfin import *
from mshr import *
import numpy as np

set_log_active(False) #Mutes the system messages

cx, cy, radius = 0.2, 0.2, 0.05
lx, ly = 2.2, 0.41

resolution = 32

class circle(SubDomain):
    def inside(self, x, on_boundary):
        return pow(x[0] - cx, 2) + pow(x[1] - cy, 2) <= pow(radius, 2) + DOLFIN_EPS

# Define 2D geometry
domain1 = Rectangle(Point(0.0, 0.0), Point(lx, ly))
domain2 = Circle(Point(cx,cy),radius, resolution)
domain = domain1 - domain2

# Generate mesh
mesh = generate_mesh(domain, resolution)

#Establish Finite Element Space
V = VectorFunctionSpace(mesh, "CG", 2, dim=2)
Q = FunctionSpace(mesh, "CG", 1)
W = V*Q
w = Function(W)
(v,q) = TestFunctions(W)
(u,p) = TrialFunctions(W)

#Problem Settings
t = 0.0 #t-start
t_max = 4.0
n_steps = 100
reynolds = 600.0
delta_t = (t_max / n_steps)
invdt = Constant((1.0 / delta_t))
nu = Constant((1 / reynolds))
fpi_count_max = 100

um = Function(V)
uold = Function(V)

zero_function = Expression(('0.0','0.0'))
inflow = Expression(('x[1]*(0.41 - x[1])*sin(pi*t*0.125)/(0.41*0.41)','0.0'),t=0.0)
zero_pressure = Expression(('0.0'))

#Define Boundaries
def right(x): return x[0] > 2.2 - DOLFIN_EPS
def left(x): return x[0] < DOLFIN_EPS
def top(x): return x[1] > 0.41 - DOLFIN_EPS
def bottom(x): return x[1] < DOLFIN_EPS
def cyl(x): return sqrt((x[0]-0.2)*(x[0]-0.2) + (x[1]-0.2)*(x[1]-0.2)) < 0.05 + DOLFIN_EPS

bc_u_left = DirichletBC(W.sub(0), inflow, left)
bc_u_top = DirichletBC(W.sub(0), zero_function, top)
bc_u_bot = DirichletBC(W.sub(0), zero_function, bottom)
bc_u_cyl = DirichletBC(W.sub(0), zero_function, cyl)
bc_u_right = DirichletBC(W.sub(1), zero_pressure, right)


bcs = [bc_u_left, bc_u_top, bc_u_bot, bc_u_cyl, bc_u_right]

NSECkNs = invdt*(inner(u,v))*dx + 0.5*nu*inner(nabla_grad(u),nabla_grad(v))*dx  - 0.5*((Dx(u[1],0)-Dx(u[0],1))*uold[1])*v[0]*dx + 0.5*((Dx(u[1],0)-Dx(u[0],1))*uold[0])*v[1]*dx + nabla_div(u)*q*dx + p*nabla_div(v)*dx

NSERHS =  invdt*(inner(um,v))*dx - 0.5*nu*inner(nabla_grad(um),nabla_grad(v))*dx + 0.5*((Dx(um[1],0)-Dx(um[0],1))*um[1])*v[0]*dx - 0.5*((Dx(um[1],0)-Dx(um[0],1))*um[0])*v[1]*dx

uold.assign(inflow)
um.assign(inflow)

for tn in range(n_steps):
    t = tn*delta_t
    inflow.t = t
    print '////// Time: ' + str(t) + '////'
    for nn in range(fpi_count_max):
        fpi_error = 1
        solve(NSECkNs == NSERHS, w, bcs)
        (u,p) = w.split(True)
        #Fixed Point Iteration Evaluation
        if(nn>0):
            fpi_vals =  u.vector().array() - uold.vector().array()
            fpi_arr = np.asanyarray(fpi_vals)
            fpi_error = abs(fpi_arr).max()
            print 'FPI Error: ' + str(fpi_error) + ' FPI Count: ' + str(nn)
        if(fpi_error < 1e-08):
            um.assign(u)
            uold.assign(u)
            break
        else:
            uold.assign(u)

plot(u)

I've used several FE packages before, and have never actually run into this issue, so I know it's user error. Any help would be greatly appreciated. Thanks!

Answer 1

How do you expect the pressure to be defined, if you do not have pressure boundary conditions?

I do not think it make sense to prescribe the velocity on the outlet. You can avoid a weak term on the outlet that way, but something very similar can be achieved by setting the pressure to zero on the outlet This is equivalent to zero normal stress.

Comments

I've actually already investigated that --it doesn't appear to change anything. I'll ammend the code above to reflect that, but it produces the same velocity plot. (Typically, I prescribe pressure at a single point somewhere on the boundary... I forgot to include it after I removed it previously.)

---

I think you are just pushing it too far with the Reynolds number. It looks OK with

Re = 60.
nsteps = 1000

---

Tracing back the citations to the orginal benchmark paper (Volker, John -- Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder) he actually forces the outflow.

I was trying to reproduce the results from a paper, and it appears they misremembered their own settings, claiming a do-nothing boundary. When I force the outflow I get precisely their results.

---

There is no such thing as a do-nothing boundary.

If you do not add a weak term for the boundary, it corresponds to zero force on that boundary - unless you prescribe a Dirichlet boundary condition, in which case the test function is zero on the boundary.

In any case the length scales for a Re=600 flow are orders of magnitude smaller than the elements in your mesh, so the results are total garbage. You will see that very clearly, when you test the convergence.

---

You're correct, but I was choosing a large length scale that produced a plot of the "issue" I was complaining about so someone wouldn't spend a half hour running the problem at the correct scale. In this case, my attempt at politeness got in my own way. Point taken. (I'm investigating an LES model, so I spend a lot of time above appropriate length scale.)

When I adjusted my settings to the ones prescribed in the Volker paper I get the results I expected to see. Thanks for your help.