CBCFLOW: Womersley velocity profile in DirichletBC

Question

Hi,

I'm trying to implement a Womersley velocity profile using cbcflow as inlet with a DirichletBC. Is this possible? The errors are inline in the code below:

# 
mu = 0.00345                           # dynamic viscosity (Ns/m2)
rho = (1.056*10**(-3))                 # blood density (g/l)
nu = mu/rho                            # kinematic viscosity

# Setup analytical solution constants
Q_womersley = 1.0

# Beta is the Poiseuille pressure drop if the flow rate is stationary Q
beta = 4.0 * nu * Q_womersley / (pi * r**4)

# Setup transient flow rate coefficients
###print "Using transient bcs."

P = 0.8  #Period
tvalues = np.linspace(0.0, P)
Qfloor, Qpeak = -0.2, 1.0
Qvalues = Q_womersley * (Qfloor + (Qpeak-Qfloor))
Q_coeffs = zip(tvalues, itertools.repeat(Qvalues,len(tvalues)))

# WOMERSLEY ANALYTICAL SOLUTION
###print "Create womersley objects..."

ua = make_womersley_bcs(Q_coeffs, mesh, 1 , nu, None, boundaries) **# cbcflow  function**

for uc in ua:
    uc.set_t(t)
    pa = Expression("-beta * x[0]", beta=1.0)
    pa.beta = beta # TODO: This is not correct unless stationary...

# Define boundary conditions
V = VectorFunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "CG", 1)
boundaries = MeshFunction("size_t", mesh, "cylinder_volume9_facet_region.xml")
u0 = Constant((0., 0., 0.))     # initial velocity

noslip  = DirichletBC(V, Constant ((0,0,0)),boundaries,3)  
#noslip = (u0,3)

## Problem: How to put ua in DirichletBC?
inflow  = DirichletBC(V, Expression(ua[0],ua[1],ua[2]), boundaries,1) 
#inflow  = (ua,1)

## TypeError: Please provide a 'str', 'tuple' or 'list' for the 'cppcode' argument.

outflow = DirichletBC(Q, Constant (0), boundaries,2) ####outflow = (pa,2)

bcu = [noslip, inflow]
bcp = [outflow]

---

If I change inflow to inflow = (ua,1) gives me an error when computing velocity because the method apply is used for DirichletBC. I don't want to change the way the problem is being solved and that's why I'm trying to define the womersley profile as a DirichletBC.

Compute tentative velocity step

        begin("Computing tentative velocity")
        b1 = assemble(L1)
        [bc.apply(A1, b1) for bc in bcu]
        solve(A1, u1.vector(), b1, "gmres", "ilu")
        end()

Pressure correction

        begin("Computing pressure correction")
        b2 = assemble(L2)
        [bc.apply(A2, b2) for bc in bcp]
        solve(A2, p1.vector(), b2, "gmres", 'amg')
        end()

Velocity correction

        begin("Computing velocity correction")
        b3 = assemble(L3)
        [bc.apply(A3, b3) for bc in bcu]
        solve(A3, u1.vector(), b3, "gmres", "ilu")
        end()

Answer 1

Hi, assuming make_womersley_bcs return a list of expressions the fix should be analogous to the following

from dolfin import *

c0 = Expression('x[0]', degree=1)
c1 = Expression('x[1]', degree=1)
cs = [c0, c1]   # This is close to what Womersley is

mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, Expression(('foo', 'bar'), foo=cs[1], bar=cs[0], degree=1), 'on_boundary')