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"),
         path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/pipe_3k.xml.gz"),
         path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/pipe_24k.xml.gz"),
         path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/pipe_203k.xml.gz"),
         path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/pipe_1611k.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.

@classmethod
def default_params(cls):
    params = NSProblem.default_params()
    params.replace(
        # Time parameters
        T=None,
        dt=1e-3,
        period=0.8,
        num_periods=1.0,
        # Physical parameters
        rho=1.0,
        mu=1.0/30.0,
        )
    params.update(
        # Spatial parameters
        refinement_level=0,
        # Analytical solution parameters
        Q=1.0,
        )
    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)
    facet_domains.set_all(3)
    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:
        uc.set_t(t)
    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:
        ucomp.set_t(t)

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 = [
        Pressure(plot_and_save),
        Velocity(plot_and_save),
        ]
    postproc = NSPostProcessor({"casedir": casedir})
    postproc.add_fields(fields)

and solves the problem

solver = NSSolver(problem, scheme, postproc)
solver.solve()