Some questions about N-S cylinder

General
1

I have studied the N-S equation cylinder example in fenics tutorial . Now i want to rewrite the example using mixed element , but the result is false (no turbulence). there is my code `

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

T = 5.0            # final time
num_steps = 5000   # number of time steps
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, 64)

# Define function spaces

elem_u = VectorElement('P', mesh.ufl_cell(), 2)
elem_p = FiniteElement('P', mesh.ufl_cell(), 2)
elem_up = MixedElement([elem_u, elem_p])
V = FunctionSpace(mesh, elem_up)

# 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'

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

# Define boundary conditions
bcu_inflow    = DirichletBC(V.sub(0), Expression(inflow_profile, degree=2), inflow)
bcu_walls     = DirichletBC(V.sub(0), Constant((0, 0)), walls)
bcu_cylinder  = DirichletBC(V.sub(0), Constant((0, 0)), cylinder)
bcp_outflow   = DirichletBC(V.sub(1), Constant(0), outflow)
bcup          = [bcu_inflow, bcu_walls, bcu_cylinder, bcp_outflow]


# Define trial and test functions
u, p = TrialFunctions(V)
v, q = TestFunctions(V)

# Define functions for solutions at previous and current time steps
up_old = Function(V)
u_old, p_old = split(up_old)
up_fem = Function(V)

# Define expressions used in variational forms
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho_f = Constant(rho)

# Define symmetric gradient
def epsilon(u):
     return sym(nabla_grad(u))

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

# Define variational problem for step 1



F_ns = rho_f*dot((u - u_old) / k, v)*dx + \
     rho_f*dot(dot(u_old, nabla_grad(u)), v)*dx + \
     inner(epsilon(u), epsilon(v))*dx -\
     p*div(v)*dx+\
     div(u)*q*dx - \
     dot(f, v)*dx 


a1 = lhs(F_ns)
L1 = rhs(F_ns)


vtkfile_u = File('n_s_cylinder/velocity.pvd')

# Time-stepping
t = 0
for n in range(num_steps):

    # Update current time
     t += dt

     solve(a1 == L1, up_fem, bcup)
     u_fem, p_fem = up_fem.split()
     # Plot solution
     if n%(num_steps/5) == 0:
          print("go, on!")
          plot(u_fem, title='Velocity')
          plt.show()
          plot(p_fem, title='Pressure')
          plt.show()
     if n == num_steps -1 :
          print("finish")
     
     vtkfile_u << (u_fem, t)

    # Update previous solution
     up_old.assign(up_fem)

Are there some examples to refer to?

2

By “turbulence” do you mean a Karman vortex street? The backward Euler / Lie splitting scheme you’ve employed is first order accurate and quite dissipative. Furthermore your element pair is not stable, consider the Taylor Hood element where the polynomials in pressure space are one degree lower than the velocity space. Personal experience is that the vortices form if you switch to a second or higher order scheme. Something like Crank Nicolson with Taylor Hood elements.

2 Likes