Incorrect oscillating flow when applying variable pressure BC (windkessel) for IPCS method

Question

Dear FEniCS community,

I have some issues obtaining a proper velocity solution, when I prescribe a pressure as boundary condition. I observe a flow, that highly oscillates around the expected flow value (numerical noise).

I have a geometry consisting of multiple outlets (Aortic Arch). At the inlet I prescribe a periodic velocity curve (based on an aortic inflow pulse). At the outlets I want a windkessel afterload, that describes the pressure-flow relation.
Something simular to the first example in: http://hplgit.github.io/fenics-mixed/doc/pub/fenics-mixed-4paper.pdf

For every timestep the new pressure is calculated from the known flow (q(n)) (fraction of inflow) and the previous pressure (p(n)).
p(n+1)_outlet_n = f(p(n), q(n))

Every timestep this value is given as Dirichlet BC for the pressure. I obtain the pressure relation that I was expecting. (see figure 1) But the flow is oscillating around the expected flow curve, such that the value between the two extrema is the value that i would expect. (see figure 2). The flow changes from direction every timestep, no matter how small I set the timestep, indicating it is numerical noise.

1
2

I expect the problem to be within the declaration of the pressure BC within the 3-steps IPCS method. I've found a somewhat simular problem here: https://fenicsproject.org/qa/11819/stokes-flow-simulation-imposing-pressure-sensible-results
But their solution (use Neuman BC instead of Dirichlet) is not an option since I am using the IPCS method.

I use the optimized IPCS solver Oasis (https://github.com/mikaem/Oasis) but the code below is a simple representation of my code in Oasis.

I hope that someone can help me with this.

Thanks in advance,

Martijn

Basic idea (within time integration loop)
(total simplified version of my code below in comments)

while t < simulationTime:

# Update velocity inlet BC
inflow_expression.t = t
inflow_expression.v_max = cur_vel
bcu_inlet    = DirichletBC(V, inflow_expression, boundaries, surf_inlet)

# Obtain previous values
pt_1 = previous_pres_of_outlet_1
pt_2 = previous_pres_of_outlet_2
qnm1_1 = previous_flow_of_outlet_1
qnm1_2 = previous_flow_of_outlet_2

# Obtain current flow
q_1 = 0.3 * flow_inlet
q_2 = 0.7 * flow_inlet

# Calculate new pressure for outlet 1 and outlet 2
pnp1_outlet_1   = ((Z_1*dt + R_1*dt + Z_1*R_1*C_1)/(R_1*C_1 + dt ))*q_1 - ((Z_1*R_1*C_1)/(R_1*C_1 + dt))*qnm1_1 + (R_1*C_1/(R_1*C_1 + dt)) * pt_1
pnp1_outlet_2   = ((Z_2*dt + R_2*dt + Z_2*R_2*C_2)/(R_2*C_2 + dt ))*q_2 - ((Z_2*R_2*C_2)/(R_2*C_2 + dt))*qnm1_2 + (R_2*C_2/(R_2*C_2 + dt)) * pt_2

# Z_1, R_1, C_1 and Z_2, R_2, C_2 are known windkessel constants

# Update pressure BC
bcp_outlet_1    = DirichletBC(Q, Constant(pnp1_1),     boundaries, surf_1)
bcp_outlet_2    = DirichletBC(Q, Constant(pnp1_2),    boundaries, surf_2)

# Solve tentative velocity step
# with no-slip and velocity inflow BC

# Solve pressure correction step
# with calculated pressure at outlets as BC

# Solve velocity update step

# update time
t += dt

Comments

Simplified version of my code:

    # Read mesh
mesh = Mesh('Aorta.xml.gz')
D = mesh.topology().dim()
boundaries = MeshFunction("size_t", mesh, D - 1, mesh.domains())
ds = ds(subdomain_data = boundaries)

# Import inflow signal
# Normally I import a aortic flow pulse from a .csv file
# which I prescribe every time step, but I use a simple 
# sinus as inflow for this example

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 2)   
Q = FunctionSpace(mesh, 'P', 1)         

# Define Trial and Test functions
u = TrialFunction(V)    
v = TestFunction(V)     
p = TrialFunction(Q)    
q = TestFunction(Q) 

# Define functions for solutions at previous and current time steps
u_n = Function(V)       
u_  = Function(V)       
p_n = Function(Q)       
p_  = Function(Q) 

# define input
v_max = 1
mu = 1
rho = 1
dt = k = 0.005
pi = 3.1415 
f = Constant((0,0,0))
nx, ny, nz = inlet_normal

# define inflow expression
inflow_expression = Expression(('nx * v_max*sin(2*pi*t)',\
                                'ny * v_max*sin(2*pi*t)',\
                                'nz * v_max*sin(2*pi*t)'),\
 degree=2, v_max=v_max, nx=nx, ny=ny, nz=nz, pi=pi, t=0)

# inlet BC
bcu_inlet    = DirichletBC(V, inflow_expression , boundaries, surf_inlet)

# No-slip BC
bcu_wall    = DirichletBC(V, Constant((0, 0, 0)), boundaries, surf_wall)

# Variatational form of tentative velocity step
F1 = \
rho*dot((u - u_n) / k, v)*dx               \
 + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
 + inner(sigma(U, p_n), epsilon(v))*dx      \
 + dot(p_n*n, v)*ds                         \
 - dot(mu*nabla_grad(U)*n, v)*ds            \
 - rho*dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Variataional form of pressure correction step 
a2 = \
 dot(nabla_grad(p), nabla_grad(q))*dx

L2 = \
 dot(nabla_grad(p_n), nabla_grad(q))*dx     \
 - (rho/k)*div(u_)*q*dx     

# Variatational form of velocity update step
a3 = \
 dot(u, v)*dx

L3 = \
 dot(u_, v)*dx                              \
 - k/rho*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

#start time loop
for i in range(int(num_steps)):

    #update time
    t += dt

    #Update inflow BC
    inflow_expression.t = t

    bcu_inlet    = DirichletBC(V, inflow_expression , boundaries, surf_inlet)

    # Define Windkessel parameters
    R_bra       = 8.150590715942891E08
    R_lcar      = 2.229839875727788E09
    R_lsub      = 1.296695638888460E09
    R_aodis     = 2.372585709880160E08

    C_bra       = 1.518294861229382E-09
    C_lcar      = 7.153876910015118E-10
    C_lsub      = 1.035902304068399E-09
    C_aodis     = 1.028382538266091E-08
    C_aodis2    = 8.00000000000000E-09

    Z_bra       = 8.232919915093829E06
    Z_lcar      = 6.709648572902069E06
    Z_lsub      = 6.516058486876683E06
    Z_aodis     = 1.556563627475424E07

    # determine pressure from previous solution at every outlet
    pt_bra = p_ * ds(surf_bra)      
    pt_bra = assemble(pt_bra) / area_bra    
    # idem for all other outlets

    # determine flow at every outlet
    q_bra = 0.1 * flow_in   # for example 10% of inflow
    # idem for all other outlets

    # Calculate pressure (3-element windkessel formulation)
    pnp1_bra   = ((Z_bra*dt + R_bra*dt + Z_bra*R_bra*C_bra)/(R_bra*C_bra + dt ))*q_bra - ((Z_bra*R_bra*C_bra)/(R_bra*C_bra + dt))*qnm1_bra + (R_bra*C_bra/(R_bra*C_bra + dt)) * pt_bra
    pnp1_lcar  = ((Z_lcar*dt + R_lcar*dt + Z_lcar*R_lcar*C_lcar)/(R_lcar*C_lcar + dt ))*q_lcar - ((Z_lcar*R_lcar*C_lcar)/(R_lcar*C_lcar + dt))*qnm1_lcar + (R_lcar*C_lcar/(R_lcar*C_lcar + dt)) * pt_lcar
    pnp1_lsub  = ((Z_lsub*dt + R_lsub*dt + Z_lsub*R_lsub*C_lsub)/(R_lsub*C_lsub + dt ))*q_lsub - ((Z_lsub*R_lsub*C_lsub)/(R_lsub*C_lsub + dt))*qnm1_lsub + (R_lsub*C_lsub/(R_lsub*C_lsub + dt)) * pt_lsub
    pnp1_aodis = ((Z_aodis*dt + R_aodis*dt + Z_aodis*R_aodis*C_aodis)/(R_aodis*C_aodis + dt ))*q_aodis - ((Z_aodis*R_aodis*C_aodis)/(R_aodis*C_aodis + dt))*qnm1_aodis + (R_aodis*C_aodis/(R_aodis*C_aodis + dt)) * pt_aodis

    # Update outlet BC
    bcp_bra     = DirichletBC(Q, Constant(pnp1_bra),     boundaries, surf_bra)
    bcp_lcar    = DirichletBC(Q, Constant(pnp1_lcar),    boundaries, surf_lcar)
    bcp_lsub    = DirichletBC(Q, Constant(pnp1_lsub),    boundaries, surf_lsub)
    bcp_aodis   = DirichletBC(Q, Constant(pnp1_aodis),   boundaries, surf_aodis)

    # make BC lists
    bcu = [bcu_inlet, bcu_wall]
    bcp = [bcp_bra, bcp_lcar, bcp_lsub, bcp_aodis]

    # Apply BC's and solve tentative velocity step
    b1 = assemble(L1)
    [bc.apply(A1, b1) for bc in bcu]
    solve(A1, u_.vector(), b1, solver, preconditioner) 

    # Apply BC's and solve pressure correction step
    b2 = assemble(L2)
    [bc.apply(A2, b2) for bc in bcp]
    solve(A2, p_.vector(), b2, solver, preconditioner)

    # solve velocity update step step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, solver)