Problems debugging sip condition for Navier Stokes


147
views
0
7 weeks ago by
Hello,

I am having a hard time to debug a simple code below in order to implement slip condition. I am trying to use slip boundary condition on the walls (north and south) for the flow around cylinder. 
I would appreciate any help.
from dolfin import *
import mshr
%matplotlib inline

# Define domain
center = Point(0.2, 0.2)
radius = 0.05
L = 2.2
W = 0.41
geometry = mshr.Rectangle(Point(0.0, 0.0), Point(L, W))  \
         - mshr.Circle(center, radius, 10)

# Build mesh
mesh = mshr.generate_mesh(geometry, 50)
plot(mesh)

# Construct facet markers
bndry = FacetFunction("size_t", mesh)
for f in facets(mesh):
    mp = f.midpoint()
    if near(mp[0], 0.0): # inflow
        bndry[f] = 1
    elif near(mp[0], L): # outflow
        bndry[f] = 2
    elif near(mp[1], 0.0) or near(mp[1], W): # walls
        bndry[f] = 3
    elif mp.distance(center) <= radius: # cylinder
        bndry[f] = 5

# Build function spaces (Taylor-Hood)
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W0 = MixedElement([V,P])
W = FunctionSpace(mesh,W0)

# No-slip boundary condition for velocity on walls and cylinder - boundary id 3
noslip = Constant((0, 0))
bc_cylinder = DirichletBC(W.sub(0), noslip, bndry, 5)

# Inflow boundary condition for velocity - boundary id 1
v_in = Expression(("0.3 * 4.0 * x[1] * (0.41 - x[1]) / ( 0.41 * 0.41 )", "0.0"),degree=2)
bc_in = DirichletBC(W.sub(0), v_in, bndry, 1)

# Collect boundary conditions
bcs = [bc_cylinder, bc_in]

# Facet normal, identity tensor and boundary measure
n = FacetNormal(mesh)
I = Identity(mesh.geometry().dim())
ds = Measure("ds", subdomain_data=bndry)
nu = Constant(0.001)

# constants :
gamma = Constant(1e2)
h     = CellSize(mesh)
eta   = Constant(1.0)
f     = Constant ((0.0,0.0))
beta  = Constant(10.0)
u_n = Constant(0.0)
def navier_stokes():

    v, q = TestFunctions(W)
    w = Function(W)
    u, p = split(w)
   
    # Define variational forms
    def sigma(u,p): return(-p*I + 2.0*nu*sym(grad(u)))
    
    t = dot(sigma(u,p),n)
    s = dot(sigma(v,q),n)
    
    B_o = + inner(sigma(u,p),grad(v))*dx - div(u)*q*dx
    B_v = + inner(grad(u)*u, v)*dx
    B_g = - dot(n,t)*dot(v,n)*ds(3) - dot(u,n)*dot(s,n)*ds(3) + gamma/h*dot(u,n)*dot(v,n)*ds(3) + beta*dot(u, v)*ds(3)
    B_w = - inner(dot(sigma(u,p),n),v)*ds(2) - inner(dot(sigma(v,q),n),u)*ds(2) + gamma/h*inner(v,u)*ds(2)

    F = dot(f,v)*dx + (gamma/h)*u_n*dot(v,n)*ds(3) - inner(dot(sigma(v,q),n),v_in)* ds(2) \
       + (gamma/h)*inner(v,v_in)*ds(2)

    solve(B_o + B_v + B_w == F, w, bcs)
    
    return w

ufile = File("stokes/velocity.pvd")
pfile = File("stokes/pressure.pvd")

if __name__ == "__main__":

    for problem in [navier_stokes]:
        begin("Running '%s'" % problem.__name__)

        # Call solver
        w = problem()

        # Extract solutions
        u, p = w.split()

        # Save solution
        u.rename("u", "velocity")
        p.rename("p", "pressure")
        ufile << u
        pfile << p

        end()

    interactive()​
Community: FEniCS Project
1
Nitsche applied to Navier equations... interesting. I wonder if your weak form has been derived correctly.
written 7 weeks ago by pf4d  
Hi pf4d!

I am almost sure the weak form is right, but I can edit and put the equations if you have time to take a look.

Thanks for your time.
written 7 weeks ago by Paulo Vicente  

1 Answer


3
7 weeks ago by
Running your code, I get
*** Error:   Unable to define linear variational problem a(u, v) == L(v) for all v.
*** Reason:  Expecting the left-hand side to be a bilinear form (not rank 1).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.​
so you're trying to define a LinearVariationalProblem for a nonlinear problem. If you change  solve(B_o + B_v + B_w == F, w, bcs) to
J = derivative(B_o + B_v + B_w - F, w)
problem = NonlinearVariationalProblem(B_o + B_v + B_w - F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()​
it is running (version 2016.2.0). Or is there any other problem that I missed?
Yes, my version is 2016.2 and you don’t miss anything. Do you know if I can use this procedure with only one iteration?

Thanks for your help!
written 7 weeks ago by Paulo Vicente  
Of course, Navier includes the acceleration term ...
written 6 weeks ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »