Based on the program to simulate the flow past a cylinder, I am currently trying to add 2 jets located at the top and the bottom of the cylinder each with a width of 10 degrees to actively control the vortex shedding, as the following picture shows.
So I need to add new boundary conditions for the jets,
jet_top = 'on_boundary && x[0]>(0.2-0.05*sin(5*pi/180)) && x[0]<(0.2+0.05*sin(5*pi/180)) && x[1]>0.2 && x[1]<0.3'
jet_bottom = 'on_boundary && x[0]>(0.2-0.05*sin(5*pi/180)) && x[0]<(0.2+0.05*sin(5*pi/180)) && x[1]<0.2 && x[1]>0.1'
and update the boundary conditions related to the jets in each time step(I haven’t completed this step, but I’ve created a function to change the boundary conditions related to the jets.
# Jet profile. Jet1 is at the top of the cylinder, and Jet2 is at the bottom of the cylinder.
jet1_f=Expression(('cos(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_1)/width)*Qjet*pi/(2*width*pow(radius,2))',\
'sin(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_1)/width)*Qjet*pi/(2*width*pow(radius,2))'),\
Qjet=0,width=Constant(10),radius=Constant(10),theta0_1=Constant(0.5*pi),degree=2)
jet2_f=Expression(('cos(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_2)/width)*(-Qjet)*pi/(2*width*pow(radius,2))',\
'sin(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_2)/width)*(-Qjet)*pi/(2*width*pow(radius,2))'),\
Qjet=0,width=Constant(10),radius=Constant(10),theta0_2=Constant(1.5*pi),degree=2)
# boundary conditions.
bcu_inflow=DirichletBC(V,inflow_f,inflow)
bcu_walls=DirichletBC(V,Constant((0,0)),walls)
bcu_cylinder=DirichletBC(V,Constant((0,0)),cylinder)
bcp_outflow=DirichletBC(Q,Constant(0),outflow)
bcp=[bcp_outflow]
bcu_jet_top=DirichletBC(V,jet1_f,jet_top)
bcu_jet_bottom=DirichletBC(V,jet2_f,jet_bottom)
bcu=[bcu_inflow, bcu_walls, bcu_cylinder,bcu_jet_top,bcu_jet_bottom]
#update the boundary conditons related to the jets
def update_jetBCs(new_Qjet):
jet1_f.Qjet=new_Qjet
jet2_f.Qjet=new_Qjet
bcu_jet_top=DirichletBC(V,jet1_f,jet_top)
bcu_jet_bottom=DirichletBC(V,jet2_f,jet_bottom)
bcu=[bcu_inflow, bcu_walls, bcu_cylinder,bcu_jet_top,bcu_jet_bottom]
return
Now I can still run this program, and the program computes the drag_coefficient and lift_coefficient as well. Yet, there are a lot of warnings “*** Warning: Found no facets matching domain for boundary condition.” Is my way of defining the new domains or boundary conditions wrong? How should I fix this problem?
Thanks for any help in advance!
The entire program is as follows. But the last version of it (without adding the lines I listed above) works well.
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import meshio
from math import sin, cos, pi
from matplotlib.animation import FuncAnimation
print("Hello1")
T=5
num_steps=5010
dt=T/num_steps
D=0.1
Re=100
U_m=1.5
mu=2*U_m*D/(3*Re)
rho=1
Q0=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)
#Refine the mesh near the cylinder
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
origin = Point(2*D, 2*D)
for cell in cells(mesh):
p = cell.midpoint()
if p.distance(origin) < 0.12:
cell_markers[cell] = True
mesh = refine(mesh, cell_markers, redistribute = True)
# 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'
jet_top = 'on_boundary && x[0]>(0.2-0.05*sin(5*pi/180)) && x[0]<(0.2+0.05*sin(5*pi/180)) && x[1]>0.2 && x[1]<0.3'
jet_bottom = 'on_boundary && x[0]>(0.2-0.05*sin(5*pi/180)) && x[0]<(0.2+0.05*sin(5*pi/180)) && x[1]<0.2 && x[1]>0.1'
# jet_top = 'on_boundary && between(x[0],[(0.2-0.05*sin(5*pi/180)),(0.2+0.05*sin(5*pi/180))],4e-16) && x[1]>0.2 && x[1]<0.3'
# jet_bottom = 'on_boundary && between(x[0],[(0.2-0.05*sin(5*pi/180)),(0.2+0.05*sin(5*pi/180))],4e-16) && x[1]<0.2 && x[1]>0.1'
# Inflow profile
inflow_profile = ('4.0*(U_m)*x[1]*(0.41-x[1])/pow(0.41,2)','0')
inflow_f=Expression(inflow_profile,U_m=Constant(1.5),degree=2)
# Jet profile. Jet1 is at the top of the cylinder, and Jet2 is at the bottom of the cylinder.
jet1_f=Expression(('cos(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_1)/width)*Qjet*pi/(2*width*pow(radius,2))',\
'sin(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_1)/width)*Qjet*pi/(2*width*pow(radius,2))'),\
Qjet=0,width=Constant(10),radius=Constant(10),theta0_1=Constant(0.5*pi),degree=2)
jet2_f=Expression(('cos(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_2)/width)*(-Qjet)*pi/(2*width*pow(radius,2))',\
'sin(atan2(x[1],x[0]))*cos(pi*(atan2(x[1],x[0])-theta0_2)/width)*(-Qjet)*pi/(2*width*pow(radius,2))'),\
Qjet=0,width=Constant(10),radius=Constant(10),theta0_2=Constant(1.5*pi),degree=2)
# boundary conditions.
bcu_inflow=DirichletBC(V,inflow_f,inflow)
bcu_walls=DirichletBC(V,Constant((0,0)),walls)
bcu_cylinder=DirichletBC(V,Constant((0,0)),cylinder)
bcp_outflow=DirichletBC(Q,Constant(0),outflow)
bcp=[bcp_outflow]
bcu_jet_top=DirichletBC(V,jet1_f,jet_top)
bcu_jet_bottom=DirichletBC(V,jet2_f,jet_bottom)
bcu=[bcu_inflow, bcu_walls, bcu_cylinder,bcu_jet_top,bcu_jet_bottom]
#update the boundary conditons related to the jets
def update_jetBCs(new_Qjet):
jet1_f.Qjet=new_Qjet
jet2_f.Qjet=new_Qjet
bcu_jet_top=DirichletBC(V,jet1_f,jet_top)
bcu_jet_bottom=DirichletBC(V,jet2_f,jet_bottom)
bcu=[bcu_inflow, bcu_walls, bcu_cylinder,bcu_jet_top,bcu_jet_bottom]
return
#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)
stress_tensor=sigma(u,p)
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_data=boundary_parts,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 = abs(2 * drag_force / (rho * 1.0 * D))
lift_coefficient = abs(2 * lift_force / (rho * 1.0 * D))
return drag_coefficient, lift_coefficient
# def compute_drag_lift_coefficients(u,p):
# n = FacetNormal(mesh) # Normal pointing out of obstacle
# 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
# Cy=MeshFunction("size_t",mesh,mesh.topology().dim()-1)
# Cy.set_all(0)
# Gamma_1=CylinderBoundary()
# Gamma_1.mark(boundary_parts,1)
# dObs = Measure("ds", domain=mesh, subdomain_id=1)
# u_t = inner(as_vector((n[1], -n[0])), u)
# drag = assemble(2 / 0.1 * (mu / rho * inner(grad(u_t), n) * n[1] - p_ * n[0]) * dObs)
# lift = assemble(-2 / 0.1 * (mu / rho * inner(grad(u_t), n) * n[0] + p_ * n[1]) * dObs)
# return drag, lift
# 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
# Create progress bar
# progress=Progress('Time-stepping')
# set_log_level(PROGRESS)
#Time-stepping
t=0
Cd_max=0
Cl_max=0
x=range(50,5010,50)
y=np.ones(100)
z=np.ones(100)
print("Hello")
pressure_data=[]
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')
drag_coefficient, lift_coefficient = compute_drag_lift_coefficients(u_,p_)
if drag_coefficient > Cd_max and n>49:
Cd_max = drag_coefficient
if lift_coefficient > Cl_max and n>49:
Cl_max = lift_coefficient
if n % 50 ==0 and n >1 :
i = n/50
y[int(i)-1]= drag_coefficient
z[int(i)-1]= lift_coefficient
print("Time step:",n,"Cd:",drag_coefficient,"Cl",lift_coefficient)
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_)
print("Cd_max = ", Cd_max, "Cl_max = ", Cl_max)
plt.plot(x,y)
plt.show()