Womersley flow in 3D

In this demo it is demonstrated how to handle problems with time-dependent boundary conditions and known analytical solution/reference solution. The problem is transient Womersley flow in a cylindrical pipe.

The source code can be found in Womersley3D.py.

We start by importing cbcflow and dolfin:

from cbcflow import *
from dolfin import *

Specifying the domain

Our domain is a cylinder of length 10.0 and radius 0.5:

LENGTH = 10.0
RADIUS = 0.5

The meshes for this has been pregenerated and is available in the demo data, see Demos.

files = [path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/pipe_1k.xml.gz"),

We define SubDomain classes for inflow and outflow:

class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < 1e-6 and on_boundary

class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > LENGTH-1e-6 and on_boundary

We could also add a SubDomain-class for the remaining wall, but this will be handled later.

Defining the NSProblem

We first define a problem class inheriting from NSProblem:

class Womersley3D(NSProblem):

The parameters of the problem are defined to give a Reynolds number of about 30 and Womersley number of about 60.

def default_params(cls):
    params = NSProblem.default_params()
        # Time parameters
        # Physical parameters
        # Spatial parameters
        # Analytical solution parameters
    return params

In the constructor, we load the mesh from file and mark the boundary domains relating to inflow, outflow and wall in a FacetFunction:

def __init__(self, params=None):
    NSProblem.__init__(self, params)

    # Load mesh
    mesh = Mesh(files[self.params.refinement_level])

    # We know that the mesh contains markers with these id values
    self.wall_boundary_id = 0
    self.left_boundary_id = 1
    self.right_boundary_id = 2

    facet_domains = FacetFunction("size_t", mesh)
    DomainBoundary().mark(facet_domains, self.wall_boundary_id)
    Inflow().mark(facet_domains, self.left_boundary_id)
    Outflow().mark(facet_domains, self.right_boundary_id)

We then define a transient profile for the flow rate, for use later:

 # Setup analytical solution constants
Q = self.params.Q
self.nu = self.params.mu / self.params.rho

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

# Setup transient flow rate coefficients
print "Using transient bcs."
P = self.params.period
tvalues = np.linspace(0.0, P)
Qfloor, Qpeak = -0.2, 1.0
Qvalues = Q * (Qfloor + (Qpeak-Qfloor)*np.sin(pi*((P-tvalues)/P)**2)**2)
self.Q_coeffs = zip(tvalues, Qvalues)

Finally, we store the mesh and facet domains to self:

# Store mesh and markers
self.initialize_geometry(mesh, facet_domains=facet_domains)

The analytical solution

The Womersley profile can be obtained by using the helper function make_womersley_bcs(). This function returns a list of scalar Expression instances defining the Womersley profile:

def analytical_solution(self, spaces, t):
    # Create womersley objects
    ua = make_womersley_bcs(self.Q_coeffs, self.mesh, self.left_boundary_id, self.nu, None, self.facet_domains)
    for uc in ua:
    pa = Expression("-beta * x[0]", beta=1.0)
    pa.beta = self.beta # TODO: This is not correct unless stationary...
    return (ua, pa)

Note that the pressure solution defined here is not correct in the transient case.

Using an analytical/reference solution

If one for example wants to validate a scheme, it is required to define the following functions:

def test_fields(self):
    return [Velocity(), Pressure()]

def test_references(self, spaces, t):
    return self.analytical_solution(spaces, t)

The test_fields() function tells that the fields Velocity and Pressure should be compared to the results from test_references(), namely the analytical solution.

These functions are used in the regression/validation test suite to check and record errors.

Initial conditions

As initial conditions we simply use the analytical solution at t=0.0:

def initial_conditions(self, spaces, controls):
    return self.analytical_solution(spaces, 0.0)

Boundary conditions

At the boundaries, we also take advantage of the analytical solution, and we set no-slip conditions at the cylinder walls:

def boundary_conditions(self, spaces, u, p, t, controls):

# Create no-slip bcs d = len(u) u0 = [Constant(0.0)] * d noslip = (u0, self.wall_boundary_id)

# Get other bcs from analytical solution functions ua, pa = self.analytical_solution(spaces, t)

# Create inflow boundary conditions for velocity inflow = (ua, self.left_boundary_id)

# Create outflow boundary conditions for pressure p_outflow = (pa, self.right_boundary_id)

# Return bcs in two lists bcu = [noslip, inflow] bcp = [p_outflow]

return (bcu, bcp)

Now, since these boundary conditions are transient, we need to use the update() function. The boundary_conditions() function is called at the start of solve step, and a call-back is done to the update() function to do any updates to for example the boundary conditions. In here, we update the time in the inlet boundary condition:

def update(self, spaces, u, p, t, timestep, bcs, observations, controls):
    bcu, bcp = bcs
    uin = bcu[1][0]
    for ucomp in uin:

Solving the problem

Finally, we initate the problem, a scheme and postprocessor

def main():
    problem = Womersley3D({"refinement_level": 2})
    scheme = IPCS_Stable()

    casedir = "results_demo_%s_%s" % (problem.shortname(), scheme.shortname())
    plot_and_save = dict(plot=True, save=True)
    fields = [
    postproc = NSPostProcessor({"casedir": casedir})

and solves the problem

solver = NSSolver(problem, scheme, postproc)