### Problems debugging sip condition for Navier Stokes

350
views
0
6 months 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)
L = 2.2
W = 0.41
geometry = mshr.Rectangle(Point(0.0, 0.0), Point(L, W))  \

# 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

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_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 6 months 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.

written 6 months ago by Paulo Vicente

3
6 months ago by
*** 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?