### Problems debugging sip condition for Navier Stokes

147

views

0

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.

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 Answer

3

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!

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.

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.