Hello there,
I am trying to implement a time-dependent incompressible Navier-Stokes flow, and it ‘sorta works’, but it gets unstable at the outlet when the velocity gets high.
I have tried to run the code on a finer mesh and with smaller time steps, but that does not change anything about the instability - it still happens at the same time.
Therefore I believe something is wrong with my formulation of Naiver-Stokes. Can someone help me find my mistake? I have attached a minimal working code where the problem can be seen.
#region: Load libraries----------------------------
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div, max_value
from progress.bar import Bar
import numpy as np
import mshr
import random
import sys
#import time
#endregion----------------------------------------
#region: define square mesh ------------------------------
domain = mshr.Rectangle(Point(0.0), Point(1,1))
mesh = mshr.generate_mesh(domain, 30)
d = mesh.topology().dim()
mf = MeshFunction("size_t",mesh,d,mesh.domains())
#endregion ----------------------------------------
#region: define boundaries ------------------------
class Top_wall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 1.0))
class Bottom_wall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 0.0))
class leftWall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 0.0))
class rightWall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 1.0))
#call boundaries
top = Top_wall()
bottom = Bottom_wall()
leftwall=leftWall()
rightwall=rightWall()
#mark boundaries
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1,0)
#names for easyness
Top=1
Bottom=2
inlet=3
outlet=4
top.mark(boundaries, Top)
bottom.mark(boundaries, Bottom)
leftwall.mark(boundaries, inlet)
rightwall.mark(boundaries, outlet)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
boundary_file = File("boundaries.pvd")
boundary_file << boundaries
#endregion ----------------------------------------
#region: define parameters ------------------------
t = 0.0; dt = 0.0256; T = 3
force = Constant((0., 0.))
nu = Constant(0.01)
mu=nu
n = FacetNormal(mesh)
#endregion----------------------------------------
#region: define functionspaces and functions-------
#define function spaces
V = FunctionSpace(mesh, 'RT', 1) # velocity space, Raviart–Thomas space
Q = FunctionSpace(mesh, 'DG', 0) # pressure space, Discontinuous Lagrange space
u1 = Function(V) #for storing velocity solution
p1 = Function(Q) #for current pressure
#define mixed functionspace
P2 = VectorElement("CG", mesh.ufl_cell(), 2)
P1 = FiniteElement("CG", mesh.ufl_cell(), 1)
element = MixedElement([P2, P1])
X = FunctionSpace(mesh, element)
# #define functions
xh = Function(X)
tentative = Function(X)
u0, p0 = split(tentative)
#define trial and test functions
(u, p) = TrialFunctions(X) #for velocity, pressure
(v,q) = TestFunctions(X)
#endregion-----------------------------------------
#region: define boundary conditions----------------
pd=Constant(1)
pd_p=Constant(0.01)
bcp1 = DirichletBC(X.sub(1), pd, boundaries, inlet)
bcp2 = DirichletBC(X.sub(1), Constant(0), boundaries, outlet)
bcu1 = DirichletBC(X.sub(0), (0,0), boundaries, Top)
bcu2 = DirichletBC(X.sub(0), (0,0), boundaries, Bottom)
bc = [bcp1, bcp2, bcu1, bcu2]
#endregion---------------------------------------
#region: define weak formulation Navier-Stokes------------
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*nu*epsilon(u) - p*Identity(len(u))
U = 0.5*(u+u0)
#IPCS method, ds over inlet and outlet
F1_ns = dot((u - u0) / dt, v)*dx \
+ dot(dot(u, nabla_grad(u0)), v)*dx \
+ inner(sigma(U, p), epsilon(v))*dx \
- dot(force, v)*dx \
+ dot(nabla_grad(p), nabla_grad(q))*dx\
+ (dot(p*n, v)- dot(nu*nabla_grad(U)*n, v))*ds(inlet) \
+ (dot(p*n, v)- dot(nu*nabla_grad(U)*n, v))*ds(outlet)
a1 = lhs(F1_ns)
L1 = rhs(F1_ns)
#endregion--------------------------------------------
#region: create save files--------------------------
#get scripts name
name_of_self = sys.argv[0].split('/')[-1].split('.')[0]
#for storing solutions
directory_xdmf='results/'+name_of_self+'/'
xdmffile_velocity = XDMFFile(directory_xdmf+'velocity.xdmf')
xdmffile_pressure = XDMFFile(directory_xdmf+'pressure.xdmf')
xdmffile_temperature = XDMFFile(directory_xdmf+'temperature.xdmf')
xdmffile_velocity.parameters["flush_output"] = True
xdmffile_pressure.parameters["flush_output"] = True
xdmffile_temperature.parameters["flush_output"] = True
#endregion-----------------------------------------
#region: solve problem---------------------------
set_log_active(False)
bar = Bar('Processing', max=T/dt)
while t +dt < T + DOLFIN_EPS:
tentative.vector()[:] = xh.vector()
u0, p0 = split(tentative)
dt = min(T-t, dt)
solve(a1 == L1, xh, bc)
u1.assign(xh.sub(0, deepcopy=True))
p1.assign(xh.sub(1, deepcopy=True))
xdmffile_velocity.write(u1, t); xdmffile_pressure.write(p1, t)
#update
t += dt
bar.next()
bar.finish()
#endregion-----------------------------------------
Thank you!
