Boundary condition the initial solution of the cylinder past flow

Question

Hi all,

I'm using the FEniCS to simulate the flow past a cylinder in the 2-D channel. The geometry and the mesh are plotted in the following picture. The length is 800 and the height is 200.

The flow is driven by a constant pressure at the inlet, where the pressure = 1. The pressure at the outlet is zero. Furthermore, the non-slip boundary conditions are adopted with regards to the two parallel walls and cylinder surface.

My code is from the modification of the demo in https://fenicsproject.org/documentation/dolfin/2016.1.0/python/demo/documented/navier-stokes/python/documentation.html. Only two pieces of code are different.

Mesh generation:

from mshr import *
# Mesh generation
r = Rectangle(Point(0, 0), Point(length, height))
c = Circle(Point(0.2 * length, 0.5 * height), radius)
g = r - c
mesh = generate_mesh(g, 40)

Boundary definition:

# Define constant boundary condition
p_in = Constant(1)
p_out = Constant(0)

# Define boundary conditions
nonslip = DirichletBC(V, (0, 0), "on_boundary && \
    (x[0] > 1 + DOLFIN_EPS && x[0] < 800 - DOLFIN_EPS)")
p_inflow = DirichletBC(Q, p_in, "x[0] <= 1 + DOLFIN_EPS")
outflow = DirichletBC(Q, p_out, "x[0] >= 800 - DOLFIN_EPS")
bcu = [nonslip]
bcp = [p_inflow, outflow]

It is expected to show the well-known vortex shedding flow. However, I only find the velocity goes larger and larger as the simulation goes on. Could someone give me some advises? Thank you in advance.

=======================================================================
October 31

Notice that the though the velocity is growing larger, it is actually very small (10^{-6}). Then I found that when I assigned an initial value for u_0, which is poiseuille and of magnitude 0.1, the velocity no longer kept growing larger. Then I added an inlet boundary condition on velocity and at the same time deleted the inlet pressure boundary condition. The result is quite unphysical---the pressure, in the domain except for the outlet, is negative.

So my new question is how to define the boundary condition and initial function properly? Are both the non-zero conditions on pressure and velocity necessary? Or is there any problems with my code?

Thank you so much.

Comments

Hi,
I think there is something wrong with the boundary conditions, the coordinates of the boundaries in particular.

---

Hi Eleonora, thank you for your suggestion. I have checked the boundary conditions by printing, within every time step, the function values of one point on the cylinder surface and one point on the wall as

print "u(160, 125) = (%e, %e), u(400, 0) = (%e, %e)"  (u1(160, 125)[0], u1(160, 125)[1], u1(400, 0)[0], u1(400, 0)[1])

It turns out that the values are keep closing to zero. Note that the magnitude of growing velocity is very small (10^{-6}). Furthermore, I found the growing velocity cannot be seen when I set a initial value of function u_0 by

u0 = Function(V)
u0.assign(u_ini)

"u_ini" is the poiseuille flow of magnitude 0.1. So I figure that the phenomena that the velocity is growing is because of the two small magnitude of the initial value.

Answer 1

Hi,

I suspect what Eleonora meant, and where at least part of the problem lies is in the line

# Define boundary conditions
nonslip = DirichletBC(V, (0, 0), "on_boundary && \
 (x[0] > 1 + DOLFIN_EPS && x[0] < 800 - DOLFIN_EPS)")

I suspect that this refers to the top and bottom wall and therefore should read

# Define boundary conditions
nonslip = DirichletBC(V, (0, 0), "on_boundary && \
 (x[1] <  DOLFIN_EPS | x[1] > 200 - DOLFIN_EPS)")

Also remember that the conditions hold for being on the boundary, not inside the region, hence the change of sign in the conditions.

Doug

Comments

Hi Doug, thank you. I hope I have gotten your meaning. Do you mean that, in my code, the line

on_boundary && (x[0] > 1 + DOLFIN_EPS && x[0] < 800 - DOLFIN_EPS)

refers actually to those inside the region but not the top and bottom wall?

I'm a newbie so the thing I think about the keyword "on_boundary" might be wrong. In my opinion, only at the top and bottom wall, the inlet and outlet or the cylinder surface can "on_boundary" be true (nonzero). Thus I used the line here mean both the wall and the cylinder surface. Is my understanding of "on_boundary" right?

Regards,

Alex

---

Hi Alex,

That's an elegant way to specify the boundaries which I didn't immediately see, so apologies.

If you specify your mesh to start at (0,0) which you have done then your left wall is, I suspect, at 0. Try the following:

# Define boundary conditions
nonslip = DirichletBC(V, (0, 0), "on_boundary && \
    (x[0] > DOLFIN_EPS && x[0] < 800 - DOLFIN_EPS)")
p_inflow = DirichletBC(Q, p_in, "x[0] <= DOLFIN_EPS")
outflow = DirichletBC(Q, p_out, "x[0] >= 800 - DOLFIN_EPS")
bcu = [nonslip]
bcp = [p_inflow, outflow]

Doug

---

Hi Doug,

That's all right. I am so glad that you are kindly discuss with me. Thank you.

Actually I used the lines you typed in the first run of my code. Then I newly add the "1" to be the boundary margin. However, this change is trivial and nearly noneffective. So I don't think the specification of the boundaries is with any problem. I guess the problem should be more mathematical, like improper boundary settings may cause undesirable numerical results, or something else, maybe.

Regards,
Alex