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

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 < 1e-6 and on_boundary

class Outflow(SubDomain):
def inside(self, x, on_boundary):
return x > 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)

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", 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
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})