Description#
In a Permeation Against Vacuum (PAV) extractor, liquid LiPb is entering a EUROFER pipe with an initial tritium concentration of \(10^{18} \ \mathrm{m^{-3}}\). The outer surface of the pipe is exposed to a perfect vacuum. The tritium concentration on this surface is assumed to be zero.
The whole domain is assumed to be at \(800 \ \mathrm{K}\).
The pipe is \(4 \ \mathrm{mm}\)-thick, \(30 \ \mathrm{cm}\)-long, and has an inner radius of \(1 \ \mathrm{cm}\).
We will assume a parabolic velocity profile:
\[
v_y = v_{y,0} \ (x - r_i) \ (x + r_i)
\]
where \(v_{y,0}=40 \ \mathrm{cm/s}\) and \(r_i\) is the inner radius.
📖 Task: compute the extraction efficiency of the PAV.
The approximate efficiency \(\eta\) can be computed as:
\[
\eta = \frac{c_\mathrm{inlet} - c_\mathrm{outlet}}{c_\mathrm{inlet}}
\]
where \(c_\mathrm{inlet}\) and \(c_\mathrm{outlet}\) are the inlet and outlet average tritium concentrations, respectively.
💡Hint: use the cylindrical system of coordinates.
Meshing function#
We will give you the meshing function create_mesh that takes several geometrical arguments.
You don’t need to modify it.
from fenics import Point, CompiledSubDomain, MeshFunction, XDMFFile, SubDomain
from mshr import Rectangle, generate_mesh
inlet_id = 1
vacuum_id = 2
outlet_id = 3
axis_id = 4
id_fluid = 5
id_pipe_walls = 6
def create_mesh(length, inner_radius, pipe_thickness, refinement=300):
"""Creates a mesh for the challenge E and writes the XDMF files
to the challenge_E folder.
Args:
length (float): length of the pipe (m)
inner_radius (float): inner radius of the pipe (m)
pipe_thickness (float): pipe thickness (m)
refinement (int, optional): refinement index, the higher the index
the finer the mesh. Defaults to 300.
"""
p1 = Point(0, 0)
p2 = Point(inner_radius, length)
fluid_rectangle = Rectangle(p1, p2)
p1 = Point(inner_radius, 0)
p2 = Point(inner_radius + pipe_thickness, length)
pipe_rectangle = Rectangle(p1, p2)
domain = fluid_rectangle + pipe_rectangle
domain.set_subdomain(1, fluid_rectangle)
domain.set_subdomain(2, pipe_rectangle)
mesh = generate_mesh(domain, refinement)
# marking physical groups (volumes and surfaces)
volume_markers = MeshFunction("size_t", mesh, mesh.topology().dim())
volume_markers.set_all(1)
tol = 1e-14
id_fluid = 5
id_pipe_walls = 6
class Fluid(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= inner_radius + tol
class Pipe(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= inner_radius - tol
fluid = Fluid()
pipe = Pipe()
# marking volumes
fluid.mark(volume_markers, id_fluid)
pipe.mark(volume_markers, id_pipe_walls)
tol = 1e-14
inlet_surface = CompiledSubDomain(
"on_boundary && near(x[1], 0, tol) && x[0] < inner_radius + tol",
tol=tol,
inner_radius=inner_radius,
)
outlet_surface = CompiledSubDomain(
"on_boundary && near(x[1], outlet_x, tol) && x[0] < inner_radius + tol",
tol=tol,
inner_radius=inner_radius,
outlet_x=length,
)
bottom_surface = CompiledSubDomain("on_boundary && near(x[0], 0, tol)", tol=tol)
top_surface = CompiledSubDomain(
"on_boundary && near(x[0], top_y, tol)",
tol=tol,
top_y=inner_radius + pipe_thickness,
)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
surface_markers.set_all(0)
inlet_surface.mark(surface_markers, inlet_id)
outlet_surface.mark(surface_markers, outlet_id)
top_surface.mark(surface_markers, vacuum_id)
bottom_surface.mark(surface_markers, axis_id)
output_file = XDMFFile("challenge_E/surface_markers.xdmf")
output_file.write(surface_markers)
output_file2 = XDMFFile("challenge_E/volume_markers.xdmf")
output_file2.write(volume_markers)
Custom AverageSurfaceCylindrical class#
For this task we will need a custom derived quantity class to compute the average value on a surface in cylindrical coordinates. For more info on custom classes, see the Task Custom classes.
import festim as F
import fenics as f
class AverageSurfaceCylindrical(F.AverageSurface):
"""
Computes the average value of a field on a given surface
int(f ds) / int (1 * ds)
ds is the surface measure in cylindrical coordinates.
ds = r dr dtheta
Args:
field (str, int): the field ("solute", 0, 1, "T", "retention")
surface (int): the surface id
"""
def __init__(self, field, surface) -> None:
super().__init__(field=field, surface=surface)
self.r = None
@property
def allowed_meshes(self):
return ["cylindrical"]
def compute(self):
if self.r is None:
mesh = (
self.function.function_space().mesh()
) # get the mesh from the function
rthetaz = f.SpatialCoordinate(mesh) # get the coordinates from the mesh
self.r = rthetaz[0] # only care about r here
# dS_z = r dr dtheta , assuming axisymmetry dS_z = theta r dr
# dS_r = r dz dtheta , assuming axisymmetry dS_r = theta r dz
# in both cases the expression with self.dx is the same
avg_surf = f.assemble(
self.function * self.r * self.ds(self.surface)
) / f.assemble(1 * self.r * self.ds(self.surface))
return avg_surf
Your answer#
We’ll start by creating the mesh.
temperature = 800 # K
length = 0.3 # m
inner_radius = 1e-2 # m
pipe_thickness = 4e-3 # m
velocity = 40e-2 # m/s
create_mesh(length, inner_radius, pipe_thickness)
We will also use the following material properties:
import h_transport_materials as htm
D_pbli = (
htm.diffusivities.filter(material="lipb")
.filter(isotope="h")
.filter(author="reiter")
)[0]
S_pbli = (
htm.solubilities.filter(material="lipb").filter(isotope="h").filter(author="aiello")
)[0]
D_eurofer = htm.diffusivities.filter(material="eurofer_97").filter(author="chen")[0]
S_eurofer = htm.solubilities.filter(material="eurofer_97").filter(author="chen")[0]
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[4], line 1
----> 1 import h_transport_materials as htm
2
3 D_pbli = (
4 htm.diffusivities.filter(material="lipb")
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/h_transport_materials/__init__.py:25
22 Rg = 8.314 * ureg.Pa * ureg.m**3 * ureg.mol**-1 * ureg.K**-1
23 avogadro_nb = 6.022e23 * ureg.particle * ureg.mol**-1
---> 25 from pybtex.database import parse_file
26 from pathlib import Path
28 bib_database = parse_file(str(Path(__file__).parent) + "/references.bib")
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/pybtex/database/__init__.py:44
42 from pybtex.errors import report_error
43 from pybtex.py3compat import fix_unicode_literals_in_doctest, python_2_unicode_compatible
---> 44 from pybtex.plugin import find_plugin
47 # for python2 compatibility
48 def indent(text, prefix):
File ~/checkouts/readthedocs.org/user_builds/festim-workshop/conda/festim1/lib/python3.11/site-packages/pybtex/plugin/__init__.py:26
2 # Copyright (c) 2006-2021 Andrey Golovizin
3 # Copyright (c) 2014 Matthias C. M. Troffaes
4 #
(...) 21 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
22 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
25 import os.path # splitext
---> 26 import pkg_resources
28 from pybtex.exceptions import PybtexError
31 class Plugin(object):
ModuleNotFoundError: No module named 'pkg_resources'
my_model = F.Simulation()
# --------------
# YOUR CODE HERE
# --------------
my_model.initialise()
We need to add the advection field. This is almost identical to the task Advection-diffusion problem but adapted to work in a multi-material case.
⚠️Do not modify the following code block.
mesh_sub = f.SubMesh(my_model.mesh.mesh, my_model.mesh.volume_markers, id_fluid)
functionspace = f.VectorFunctionSpace(mesh_sub, "CG", 1)
velocity_expr = f.Expression(
("0", "v*(x[0] - inner_radius)*(x[0] + inner_radius)"),
v=velocity,
inner_radius=inner_radius,
degree=2,
)
velocity_function = f.interpolate(velocity_expr, functionspace)
V = f.VectorFunctionSpace(my_model.mesh.mesh, "CG", 1)
u = f.Function(V)
v = f.TestFunction(V)
form = f.inner(u, v) * my_model.mesh.dx
form += f.inner(velocity_function, v) * my_model.mesh.dx(id_fluid)
f.solve(form == 0, u, bcs=[])
velocity_function = u
XDMFFile("challenge_E/velocity_field.xdmf").write(velocity_function)
hydrogen_concentration = my_model.h_transport_problem.mobile.mobile_concentration()
test_function_mobile = my_model.h_transport_problem.mobile.test_function
advection_term = f.inner(
f.dot(f.grad(hydrogen_concentration), velocity_function), test_function_mobile
) * my_model.mesh.dx(id_fluid)
my_model.h_transport_problem.F += advection_term
my_model.run()
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[6], line 1
----> 1 mesh_sub = f.SubMesh(my_model.mesh.mesh, my_model.mesh.volume_markers, id_fluid)
3 functionspace = f.VectorFunctionSpace(mesh_sub, "CG", 1)
5 velocity_expr = f.Expression(
6 ("0", "v*(x[0] - inner_radius)*(x[0] + inner_radius)"),
7 v=velocity,
8 inner_radius=inner_radius,
9 degree=2,
10 )
AttributeError: 'NoneType' object has no attribute 'mesh'
# YOUR CODE GOES HERE
print(f"Efficiency: ")
Efficiency: