Hi there!
I have been working around with the formulation of Navier-Stokes using Chorin’s splitting method:
I know the code is old, but with a few fixes it runs fine. I am however wondering about some instabilities: The velocity keeps increasing, even if I set the pressure to be constant 1.
Is this because we in reality need Neuman conditions on the momentum equation, or is it something else?
Thanks.
As you have not supplied your code, there is little I can do in terms of analyzing the resulting problem.
Note that you can find the updated version here: https://bitbucket.org/fenics-project/dolfin/src/master/python/demo/documented/navier-stokes/demo_navier-stokes.py
I experience the same problem with the updated version you just send.
Here is my code:
######---------------------------------------------######
#Trying to fix wrong velocity for Navier-Stokes#
#Mesh is box#
######---------------------------------------------######
#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
from tqdm import tqdm
import sys as system
#endregion----------------------------------------
#region: Choose script conditions------------------
#get rid of solving messages
set_log_active(False)
#set mesh
generate_mesh=True
#get scripts name
name_of_self = system.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_saturation = XDMFFile(directory_xdmf+'saturation.xdmf')
#endregion-----------------------------------------
#region: Set parameters----------------------------
#set time
t=0.0; T = 10; dt = 0.0256
#parameters
k = Constant(dt)
force = Constant((0., 0.))
nu = Constant(0.01)
#endregion ----------------------------------------
#region: define mesh ------------------------------
# Create mesh
mesh = UnitSquareMesh(64, 64)
n = FacetNormal(mesh)
#endregion ----------------------------------------
#region: define boundaries
class Top(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 1.0))
class Bottom(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))
top=Top()
bottom = Bottom()
leftwall=leftWall()
rightwall=rightWall()
boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1,0)
top.mark(boundaries, 1)
bottom.mark(boundaries,2)
leftwall.mark(boundaries, 3)
rightwall.mark(boundaries, 3)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
#endregion
#region: define functions -------------------------
#define vector-vaued functionspace and scalar-valued function space
V = VectorFunctionSpace(mesh, "CG", 2) #for velocity
Q = FunctionSpace(mesh, "CG", 1) #for pressure
#define test and trial functions
v = TestFunction(V) #for velocity
u = TrialFunction(V) #for velocity
q = TestFunction(Q) #for pressure
p = TrialFunction(Q) #for pressure
# Create functions
u0 = Function(V) #for tentative velocity
u1 = Function(V) #for storing velocity solution
p1 = Function(Q) #for pressure
#endregion ----------------------------------------
#region: define boundary conditions ---------------
pd=Constant(1)
#pd=Expression('10/(1+exp(-8*t+5))',t=0.0,degree=1)
noslip = DirichletBC(V, (0, 0), boundaries,3)
inflow = DirichletBC(Q, pd, boundaries,1)
outflow = DirichletBC(Q, 0, boundaries,2)
bcu = [noslip]
bcp = [inflow, outflow]
#endregion ----------------------------------------
#region: define weak formulation --
# Tentative velocity step
F1 = (1/k)*inner(u - u0, v)*dx \
+ inner(grad(u0)*u, v)*dx \
+ nu*inner(grad(u), grad(v))*dx - inner(force, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)
# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*inner(div(u1),q)*dx
# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
# Assemble matrices
A1 = assemble(a1); b1 = None
A2 = assemble(a2); b2 = None
A3 = assemble(a3); b3 = None
#endregion ----------------------------------------
#region: solve problem ----------------------------
#time step and update
bar = Bar('Processing', max=T/dt)
while t < T + DOLFIN_EPS:
pd.t = t+dt
dt = min(T-t, dt)
#update
bar.next()
# Compute tentative velocity step
b1 = assemble(L1, tensor=b1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres")
# Pressure correction
b2 = assemble(L2, tensor=b2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2, "gmres")
# Velocity correction
b3 = assemble(L3, tensor=b3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres")
#save solutions
xdmffile_velocity.write(u1, t); xdmffile_pressure.write(p1, t)
t += dt
# Update previous solution
u0.assign(u1)
bar.finish()
#endregion ----------------------------------------I have posted my code now. Can you tell me with my issue is, thank you? (:
You are trying to find the velocity field in a rectangular channel driven by a pressure gradient and you are trying to find a steady solution to your problem. That means you need to compare your solution to the solution of the steady Poiseuille flow equation in a rectangular channel: d^2u/dx^2 = 1/mu*dp/dy. With no-slip boundary condition the solution is:
u(x) = 1/(2*mu)*dp/dy*x*(x-h)
you set p = 1 on inlet and p = 0 on outlet and use the unitsquare mesh, thus dp/dy = 1, and also h = 1. You use nu = 0.01 and rho = 1 (by omitting rho), therefore mu = nu*rho = 0.01. Inserting these numbers, the solution becomes.
u(x) = 50*x*(x-h)
Max_x | u(x) | = 50*0.5*0.5 = 12.5
Therefore you would expect the maximal velocity in the channel to be 12.5. In your example it has not even reached 8 so there is no problem with what you are describing or posting. If I run your code with dt = 1 and T = 200 I reach a steady-state maximum value of 13.5. At this point the velocity is constant. I guess the solution is a bit inaccurate due to the mesh, the coarse time step or similar. Let me know if this helps.
Thank you, I reached a similar answer.