# Steady Stokes flow converge to a false time step length dependent value using fractional step methd

Hello everyone,

I was trying to solve a steady pressure driven Stokes channel flow by using the Chorin scheme with neglecting the innertial term. But surprisingly I found the velocity field converge to a false time step length dependent value. I list my MWE below, with which I found when dt=0.1 , 0.2, 0.4, u_max=1.1, 1.2, 1.4 respectively. Did i miss some fundamental modification that required for solving a steady problem with false time approaching? Thanks!

``````from dolfin import *

# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;

N = 32
mesh = UnitSquareMesh(N, N)

# Inflow boundary
class InflowBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS

# Outflow boundary
class OutflowBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1 - DOLFIN_EPS

# No-slip boundary
class NoslipBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# Define function spaces (P2-P1), but not mixed space
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)

# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)

# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)
# ------------------------------------------------------------------------
# Define boundary conditions
# ------------------------------------------------------------------------

# Define pressure boundary condition
def pressure_bc(Q):
element = FiniteElement("CG", triangle, 1)
return Expression("1 - x[0]", element=element)

# no-slip
bv = DirichletBC(V, Constant((0.0, 0.0)), NoslipBoundary())

# Create boundary conditions for pressure
bp0 = DirichletBC(Q, pressure_bc(Q), InflowBoundary())
bp1 = DirichletBC(Q, pressure_bc(Q),  OutflowBoundary())

bcu   = [bv]
bcp   = [bp0, bp1]

# Set parameter values
dt = 0.4
T = 30
#nu = 0.01
nu = Constant(1.0/8.0)

# Define coefficients
k = Constant(dt)
f = Constant((0, 0))

# Tentative velocity step
#F1 = (1/k)*inner(u - u0, v)*dx + \
#a1 = lhs(F1)
#L1 = rhs(F1)
L1 = (1/k)*inner(u0, v)*dx + inner(f, v)*dx

# Pressure update
L2 = -(1/k)*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)
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

# Time-stepping
#t = dt
t = 0.0
while t < T + DOLFIN_EPS:

# Update pressure boundary condition
#p_in.t = t

# Compute tentative velocity step
print("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "default")

# Pressure correction
print("Computing pressure correction")
b2 = assemble(L2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2, "cg", prec)

# Velocity correction
print("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "default")

# calculate the norm of difference between two time step
u_delta = u1.vector() - u0.vector()
print("error nomr =", norm(u_delta, 'linf'))

# Move to next time step
u0.assign(u1)
t += dt
print("t =", t)

print("u_max", u1.vector().max())``````

In the first iteration you should have t=dt instead of t=0.00.
Also did you try getting rid of the Boundary conditions in the second step? These will usually cause a higher error in the helmholtz-hodge-decomposition which could effect the velocity due to step 3.

My answer is late, but maybe it still helps.