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()