Obtaining Incorrect Results for Flow Over a Cylinder - Navier Stokes

variational formulation
1

Hello! I am currently solving the incompressible flow over a cylinder problem based on the tutorial reference Solving PDEs in Python - <br> The FEniCS Tutorial Volume I . I changed some of ways this problem is solved to add additional output results (this is shown in my code below). I am verifying my results with the benchmark site DFG benchmark 2D-2 (RE100, periodic) - Featflow which simulates the same problem. I am wondering if the FEniCS tutorial is not accounting for the boundary conditions in the way this above problem is classifying them. If anyone can help me understand if I am solving this problem correctly based on the benchmarks conditions it would be greatly appreciated.

# Constant Variables
T = 5.0            # final time - Change This Number If Needed
num_steps = 5000   # number of time steps - Change This Number If Needed
dt = T / num_steps # time step size
mu = 0.001         # dynamic viscosity
rho = 1            # density

# Create Mesh
channel = Rectangle(Point(0, 0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2, 0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain, 75) # For Mesh Refinement - Change This Number

# Function Space
V = VectorFunctionSpace(mesh, 'CG', 2)
Q = FunctionSpace(mesh, 'CG', 1)

# Functions
def epsilon(u):
    return sym(nabla_grad(u))
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(len(u))

# Boundary Conditions
class inflow(SubDomain):
    def inside(self, x, on_boundary):
	    return near(x[0], 0) and on_boundary
class outflow(SubDomain):
    def inside(self, x, on_boundary):
	    return near(x[0], 2.2) and on_boundary
class walls(SubDomain):
    def inside(self, x, on_boundary):
	    return (near(x[1], 0) or near(x[1], 0.41)) and on_boundary
class cylinder(SubDomain):
    def inside(self, x, on_boundary):
	    return x[0]>=0.15 and x[0]<=0.25 and x[1]>=0.15 and x[1]<=0.25 and on_boundary
	
inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')

bc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
bc.set_all(0)
inflow().mark(bc, 1)
outflow().mark(bc, 2)
walls().mark(bc, 3)
cylinder().mark(bc, 4)

bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow()) 
bcp_outflow = DirichletBC(Q, Constant(0), outflow())
bcu_walls = DirichletBC(V, Constant((0, 0)), walls())
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder())

bcu = [bcu_inflow, bcu_walls, bcu_cylinder]
bcp = [bcp_outflow]
ds = fe.ds(subdomain_data=bc)
    
# Trial/Test Functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Time-Step Function Solutions
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Variational Form
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
I = Identity(mesh.geometry().dim())
f  = Constant((0, 0))
k  = Constant(dt)
rho = Constant(rho)
mu = Constant(mu)
nu = Constant(0.001)

# Step 1
F1 = rho*dot((u - u_n) / k, v)*dx + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx + inner(sigma(U, p_n), epsilon(v))*dx + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble Meshes
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply Boundary Conditions to Meshes
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for Visualization Output
xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')
xdmffile_u.parameters["flush_output"] = True
xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')
xdmffile_p.parameters["flush_output"] = True

# Create progress bar
progress = Progress('Time-stepping', num_steps)

# Solve Problem w/ Time-stepping
t = 0
for j in range(num_steps):

    # Update current time
    t += dt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcu]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, 'cg', 'sor')

    # Save Solution to File (XDMF/HDF5)
    xdmffile_u.write(u_, t)
    xdmffile_p.write(p_, t)

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

    # Update progress bar
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)

    a_1 = Point(0.15, 0.2)
    a_2 = Point(0.25, 0.2)
    p_diff = p_(a_1) - p_(a_2)

    I = Identity(u_.geometric_dimension())
    force = dot(-p_*I + 2.0*nu*sym(grad(u_)), n)
    D = (-20*force[0])*ds(4)
    L = (20*force[1])*ds(4)
2

A similar benchmark problem has for instance been covered in:

Using the MultiMesh functionality in FEniCS.

I have also implemented the benchmark in the dolfin-x tutorial:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code2.html

This might give you some additional insight.
Just with a quick glance, i cannot see any particular issues with your code. Please note that if you are starting from 0 flow, and introducing a sudden inflow, the solution will be unstable for the first time steps (thus why the links above do the benchmark with a gradual increase of inflow).

2 Likes
3

@ria01 RG,
I was looking at the ft08_navier_stokes_cylinder.py file and noticed a few slight differences between it and your posted code. Sometimes these small things make or break code, so here are those differences:

ft08 has function space family arguments as ‘P’ while you are using a Continuous Galerkin (standard Lagrange family) ‘CG’. (I didn’t find what ‘P’ actually was, just noted that it was used in the original ft08)

Also ft08 has:

cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'

while you have:

The difference I am noticing is the numbers used, not the Python function/C++ string syntax swap, as you seem to have done that part correctly according to the book ~but I don’t know C++, so maybe you do want to check your function! The tutorial pointed out that the idea was to find boundary points within the 2D domain.

Finally, I guess that the additional output you mention is this:

Since it wasn’t in the tutorial for ft08, and I am not familiar enough with what you want to show, I would just caution you to double-check what you are wanting to be calculated here, as my ft08 runs (I Ieft in a print statement for umax and the time) and gets a u max=2.16 at t=5.00s…

And I am still new at this, so although the benchmark has a nice new zip file with the velocity plot, I am unsure how to access it to find out if my umax is correct given the ft08 inputs…but sometimes a noob can ask a question that leads to being able to correct your code! So, here’s hoping that the discrepancies I noticed can help you debug your situation.