Getting 0 drag coefficient and lift coefficient when working on the flow over a cylinder

General
1

I 'm currently trying to get the drag coefficient and lift coefficient of the flow over a cylinder, yet the results I get are not in line with the experimental data (Cd should be around 3.2 with Re = 100 in this case). No matter the timestep, the value of both the drag and lift coefficients are 0!

Can anybody help me point out what is wrong in my code? Should I set a finer mesh or the boundary conditions near the cylinder boundary should no be the noslip bc? And how should I modify the code to get the correct result? (Does it have sth to do with the D’Alembert Paradox and how should I fix it?)

My code is as follows. Many thanks for any generous help in advance!

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import meshio

print("Hello1")

T=5
num_steps=5001
dt=T/num_steps
D=0.1
Re=100
U_m=1.5
mu=2*U_m*D/(3*Re)
rho=1


#Create mesh
channel = Rectangle(Point(0,0),Point(22*D,4.1*D))
cylinder=Circle(Point(2*D,2*D),0.5*D)
domain=channel-cylinder
mesh=generate_mesh(domain,64)


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

#Define boundaries
inflow = 'near(x[0],0)'
outflow = 'near(x[0],2.2)'
walls = 'near(x[1],0)||near(x[1],0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1] >0.1 && x[1] <0.3'

# Inflow profile
inflow_profile = ('4.0*1.5*x[1]*(0.41-x[1])/pow(0.41,2)','0')

# BCs
bcu_inflow=DirichletBC(V,Expression(inflow_profile,degree=2),inflow)
bcu_walls=DirichletBC(V,Constant((0,0)),walls)
bcu_cylinder=DirichletBC(V,Constant((0,0)),cylinder)
bcp_outflow=DirichletBC(Q,Constant(0),outflow)
bcu=[bcu_inflow, bcu_walls, bcu_cylinder]
bcp=[bcp_outflow]

#Trial and Test functions
u=TrialFunction(V)
v=TestFunction(V)
p=TrialFunction(Q)
q=TestFunction(Q)

# Functions for solutions at previous and current time steps
u_n=Function(V)
u_=Function(V)
p_n=Function(Q)
p_=Function(Q)

# Expressions used in variational forms
U=0.5*(u_n+u)
n=FacetNormal(mesh)
f=Constant((0,0))
k=Constant(dt)
mu=Constant(mu)

# Symmetric gradient
def epsilon(u):
    return sym(nabla_grad(u))

# Stress tensor
def sigma(u,p):
    return  2*mu*epsilon(u)-p*Identity(len(u))

def compute_drag_lift_coefficients(u,p):
    # Define normal vector along the cylinder surface
    n = FacetNormal(mesh) 
    stress_tensor=sigma(u,p_n)
    
    boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundary_parts.set_all(0) 
    
    class CylinderBoundary(SubDomain):
        def inside(self,x,on_boundary):
            tol=1E-14
            return on_boundary and x[0]>0.1 and x[0]<0.3 and x[1] >0.1 and x[1] <0.3
    
    
    Gamma_1=CylinderBoundary()
    Gamma_1.mark(boundary_parts,1)
    
    ds=Measure('ds',domain=mesh, subdomain_id=1)

   
    
    force = dot(stress_tensor,n)
    drag_force=assemble(force[0]*ds)
    lift_force=assemble(force[1]*ds)
    # Compute drag and lift coefficients
    drag_coefficient = 2 * drag_force / (rho * 1.0 * D)
    lift_coefficient = 2 * lift_force / (rho * 1.0 * D)

    return drag_coefficient, lift_coefficient



# Variational problem for 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)
# Variational problem for 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
# Variational problem for step 3
a3=dot(u,v)*dx
L3=dot(u_,v)*dx-k*dot(nabla_grad(p_-p_n),v)*dx

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

#Apply bcs to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u =XDMFFile('NS_cylinder/velocity.xdmf')
xdmffile_p =XDMFFile('NS_cylinder/pressure.xdmf')

# Create time series (for use in reaction_system.py)
timeseries_u=TimeSeries('NS_cylinder/velocity_series')
timeseries_p=TimeSeries('NS_cylinder/pressure_series')

# Save mesh to file
File('NS_cylinder/cylinder.xml.gz') << mesh

#Time-stepping
t=0
Cd_max=0
Cl_max=0
x=range(0,51)
y=np.ones(51)
print("Hello")


for n in range(num_steps):
    t+=dt
     #1 Tentative velocity step
    b1=assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')
    
    #2 pressure correction step
    b2=assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg' )
    
    #3 Velocity correction step
    b3=assemble(L3)
    solve(A3, u_.vector(), b3, 'cg','sor')
    
    
    
    xdmffile_u.write(u_,t)
    xdmffile_p.write(p_,t)
    
    timeseries_u.store(u_.vector(),t)
    timeseries_p.store(p_.vector(),t)
    
    u_n.assign(u_)
    p_n.assign(p_)
    
    drag_coefficient, lift_coefficient = compute_drag_lift_coefficients(u_,p_)
    if drag_coefficient > Cd_max and n>1:
        Cd_max = drag_coefficient
    if lift_coefficient > Cl_max and n>1:
        Cl_max = lift_coefficient
    
    if n % 100 ==0 and n >1 :
        i = n/100
        y[int(i)]= drag_coefficient
        
        print("Time step:",n,"Cd:",drag_coefficient,"Cl",lift_coefficient)
        
    
print("Cd_max = ", Cd_max, "Cl_max = ", Cl_max)

plt.plot(x,y)
plt.show()
2

You are not supplying the boundary markers here,
thus print(assemble(1*ds)) returns 0. change this to

    ds=Measure('ds',domain=mesh, subdomain_data=boundary_parts, subdomain_id=1)