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

variational formulation
1

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;

# Load mesh from file
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 + \
#     nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
#a1 = lhs(F1)
#L1 = rhs(F1)
a1 = (1/k)*inner(u, v)*dx + nu*inner(grad(u), grad(v))*dx 
L1 = (1/k)*inner(u0, v)*dx + inner(f, v)*dx

# Pressure update
a2 = inner(grad(p), grad(q))*dx
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())
2

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. :slight_smile: