Steady State Transport in Fluid and Pipe Permeator Wall
General Case Description
This verification case compares the TMAP8 solution with an independent analytical implementation of the steady state permeator against vacuum (PAV) transport model based on Fuerst et al. (2023). The model represents hydrogen isotope transport from flowing PbLi through a cylindrical vanadium membrane and into an ideal vacuum.
The transport process includes liquid mass transfer, jump between PbLi and vanadium surface, radial diffusion through the vanadium membrane, and recombination at the vacuum surface. The calculation is steady state and isothermal, and the hydrogen isotope partial pressure on the vacuum side is assumed to be zero.
TMAP8 represents the 5 m permeator with 20 axial control volumes. The outlet concentration of each control volume is used as the inlet concentration of the following control volume. The same radial transport model is applied within each control volume, as illustrated in Figure 1. The independent Python reference uses the same radial transport model but solves the corresponding continuous axial mass balance equation.
Figure 1: Schematic of the 20 axial control volumes and the radial transport model applied within each control volume.
Case Setup
The geometry and operating parameters are summarized in Table 1.
Table 1: Geometry and operating parameters used for the PAV verification case.
| Parameter | Description | Value | Units |
|---|---|---|---|
| Temperature | 673.15 | K | |
| Active permeator length | 5.00 | m | |
| Inner tube radius | m | ||
| Outer tube radius | m | ||
| Number of TMAP8 axial control volumes | 20 | - | |
| Inlet dissolved-isotope concentration | 1.0 | mol/m | |
| PbLi velocity | 1.78 | m/s | |
| Reynolds number | - |
The gas constant is J/(mol K), consistent with PhysicalConstants.
Material properties and liquid mass transfer
The PbLi hydrogen isotope diffusivity is (Fukada et al., 2013)
and the PbLi density is
The effective PbLi molar concentration is
where kg/mol.
The PbLi solubility coefficient is calculated using the model reported by Fuerst et al. (2023), based on Aiello et al. (2006):
with Pa.
For the vanadium membrane, the hydrogen isotope diffusivity is based on Völkl and Alefeld (1975),
and the hydrogen isotope solubility is based on Veleckis and Edwards (1969):
The hydraulic diameter is
Using the specified Reynolds number, the kinematic viscosity and Schmidt number are
and
The Sherwood number and liquid mass transfer coefficient are calculated using the correlation of Linton and Sherwood (1950):
and
The material property model parameters used in the reference implementation are summarized in Table 2. The corresponding material properties and transport quantities evaluated at 673.15 K are summarized in Table 3.
Table 2: Material-property parameters for the PAV verification case.
| Parameter | Description | Value | Units | Basis |
|---|---|---|---|---|
| PbLi diffusivity pre-exponential factor | m/s | Fukada et al. (2013) | ||
| PbLi diffusivity activation energy | J/mol | Fukada et al. (2013) | ||
| PbLi density constant | 10520.35 | kg/m | Fuerst et al. (2023) | |
| PbLi density linear coefficient | 1.19051 | kg/(m K) | Fuerst et al. (2023) | |
| PbLi solubility coefficient factor | Pa | Fuerst et al. (2023) | ||
| Vanadium diffusivity pre-exponential factor | m/s | Völkl and Alefeld (1975) | ||
| Vanadium diffusivity activation energy | J/mol | Völkl and Alefeld (1975) | ||
| Vanadium solubility pre-exponential factor | 0.138 | mol/(m Pa) | Veleckis and Edwards (1969) | |
| Vanadium solubility exponential coefficient | J/mol | Veleckis and Edwards (1969) | ||
| Vanadium surface recombination coefficient | m/(mol s) | Fuerst et al. (2021) | ||
| Sherwood correlation coefficient | 0.023 | - | Linton and Sherwood (1950) | |
| Reynolds number exponent | 0.83 | - | Linton and Sherwood (1950) | |
| Schmidt number exponent | - | Linton and Sherwood (1950) | ||
| Schmidt number | 76.0233 | - | Calculated | |
| Sherwood number | 1376.29 | - | Calculated |
Table 3: Material properties and transport quantities evaluated at 673.15 K for the PAV verification case.
| Parameter | Description | Value | Units | Basis |
|---|---|---|---|---|
| Isotope diffusivity in PbLi | m/s | Calculated | ||
| PbLi density | kg/m | Calculated | ||
| PbLi solubility coefficient | mol/(m Pa) | Calculated | ||
| Isotope diffusivity in vanadium | m/s | Calculated | ||
| Vanadium solubility coefficient | mol/(m Pa) | Calculated | ||
| Liquid mass transfer coefficient | m/s | Calculated |
Analytical solution
The analytical solution is based on the analysis by Fuerst et al. (2023). For a local bulk PbLi concentration , the liquid mass transfer relation is
where is the PbLi concentration adjacent to the membrane. Equilibrium partitioning gives
Defining
radial diffusion through the cylindrical vanadium membrane is
The membrane resistance is therefore
At the vacuum surface, the model uses
Defining
and eliminating the intermediate interface concentrations gives
The positive root is
and the local permeation flux is
Continuous axial reference solution
TMAP8 applies the radial permeation relation in 20 axial control volumes. For the analytical calculation, the axial balance is instead treated continuously:
where
Let
From the radial solution,
Combining this relation with the continuous axial balance and integrating from the inlet to axial location gives
where
To compare TMAP8 predictions to the analytical solution, we solve this scalar analytical relation at the same axial locations used by the TMAP8 control volume model and recovers the bulk concentration from
The extraction efficiency is calculated as
Limiting case of negligible surface recombination resistance
The present formulation includes finite recombination kinetics at the vacuum surface through . The analytical model of Humrickhouse and Merrill (2015) instead assumes that the dissolved isotope concentration at the outer membrane surface is approximately zero because of the vacuum, . Within the present recombination boundary condition, this assumption is recovered in the limit
for which while the permeation flux remains finite. The membrane transport relation then reduces to
Combining this expression with the liquid mass transfer relation gives
where
Using , , and the solid permeability , this becomes
which is the dimensionless radial transport parameter defined by Humrickhouse and Merrill. Substituting the limiting flux into the continuous axial balance,
with gives
Integrating over the active permeator length yields
where
Therefore, in the limit of negligible surface recombination resistance,
which is the permeator efficiency expression derived by Humrickhouse and Merrill. Thus, the finite recombination model used here contains their ideal-vacuum analytical solution as the limiting case.
Results
Figure 2 compares the 20 segment TMAP8 solution with the continuous analytical reference. The concentration profile root mean square percentage error (RMSPE) is approximately 0.127 %, which signals good agreement. As shown in Figure 2, TMAP8 starts to slightly underestimate the analytical solution at higher axial positions, which can be alleviated by increasing the number of axial control volumes and therefore better approximate the continuous solution.
Figure 2: Comparison of the 20 segment TMAP8 bulk concentration profile with the continuous analytical PAV reference solution.
The numerical comparison is summarized in Table 4, which shows acceptable performance for this verification.
Table 4: Comparison of the TMAP8 solution with the analytical reference.
| Quantity | TMAP8 | Analytical | Relative error |
|---|---|---|---|
| Outlet concentration, mol/m | 0.759275 | 0.760904 | 0.214 % |
| Extraction efficiency | 0.240725 | 0.239096 | 0.681 % |
| Average permeation flux, mol/(m s) | 0.681 % |
Input files
The input file for this case can be found at (test/tests/ver-1p/ver-1p.i), which is used as a test in TMAP8 at (test/tests/ver-1p/tests). The independent analytical comparison and figure generation are implemented in (test/tests/ver-1p/comparison_ver-1p.py).
References
- A. Aiello, A. Ciampichetti, and G. Benamati.
Determination of hydrogen solubility in lead lithium using sole device.
Fusion Engineering and Design, 81(1):639–644, 2006.
Proceedings of the Seventh International Symposium on Fusion Nuclear Technology.
URL: https://doi.org/10.1016/j.fusengdes.2005.06.364, doi:10.1016/j.fusengdes.2005.06.364.[Export]
BibTeX
@article{AIELLO2006639, author = "Aiello, A. and Ciampichetti, A. and Benamati, G.", title = "Determination of hydrogen solubility in lead lithium using sole device", journal = "Fusion Engineering and Design", volume = "81", number = "1", pages = "639--644", year = "2006", note = "Proceedings of the Seventh International Symposium on Fusion Nuclear Technology", issn = "0920-3796", doi = "10.1016/j.fusengdes.2005.06.364", url = "https://doi.org/10.1016/j.fusengdes.2005.06.364", keywords = "Lead lithium, Solubility" }RIS
TY - JOUR AU - Aiello, A. AU - Ciampichetti, A. AU - Benamati, G. TI - Determination of hydrogen solubility in lead lithium using sole device JO - Fusion Engineering and Design PY - 2006 VL - 81 IS - 1 DO - 10.1016/j.fusengdes.2005.06.364 UR - https://doi.org/10.1016/j.fusengdes.2005.06.364 SP - 639 EP - 644 ER -Plain Text
A. Aiello, A. Ciampichetti, and G. Benamati. Determination of hydrogen solubility in lead lithium using sole device. Fusion Engineering and Design, 81(1):639–644, 2006. Proceedings of the Seventh International Symposium on Fusion Nuclear Technology. URL: https://doi.org/10.1016/j.fusengdes.2005.06.364, doi:10.1016/j.fusengdes.2005.06.364. - Thomas F. Fuerst, Matthew D. Eklund, John A. Leland, Adriaan A. Riet, and Chase N. Taylor.
Parametric Study of the Vacuum Permeator for the Tritium Extraction eXperiment.
Fusion Science and Technology, 79(8):1224–1234, 2023.
URL: https://doi.org/10.1080/15361055.2023.2196237, doi:10.1080/15361055.2023.2196237.[Export]
BibTeX
@article{fuerst2023parametric, author = "Fuerst, Thomas F. and Eklund, Matthew D. and Leland, John A. and Riet, Adriaan A. and Taylor, Chase N.", title = "{Parametric Study of the Vacuum Permeator for the Tritium Extraction eXperiment}", journal = "Fusion Science and Technology", volume = "79", number = "8", pages = "1224--1234", year = "2023", publisher = "Taylor \\& Francis", doi = "10.1080/15361055.2023.2196237", url = "https://doi.org/10.1080/15361055.2023.2196237" }RIS
TY - JOUR AU - Fuerst, Thomas F. AU - Eklund, Matthew D. AU - Leland, John A. AU - Riet, Adriaan A. AU - Taylor, Chase N. TI - Parametric Study of the Vacuum Permeator for the Tritium Extraction eXperiment JO - Fusion Science and Technology PY - 2023 VL - 79 IS - 8 PB - Taylor & Francis DO - 10.1080/15361055.2023.2196237 UR - https://doi.org/10.1080/15361055.2023.2196237 SP - 1224 EP - 1234 ER -Plain Text
Thomas F. Fuerst, Matthew D. Eklund, John A. Leland, Adriaan A. Riet, and Chase N. Taylor. Parametric Study of the Vacuum Permeator for the Tritium Extraction eXperiment. Fusion Science and Technology, 79(8):1224–1234, 2023. URL: https://doi.org/10.1080/15361055.2023.2196237, doi:10.1080/15361055.2023.2196237. - Thomas F. Fuerst, Paul W. Humrickhouse, Chase N. Taylor, and Masashi Shimada.
Surface effects on deuterium permeation through vanadium membranes.
Journal of Membrane Science, 620:118949, 2021.
URL: https://doi.org/10.1016/j.memsci.2020.118949, doi:10.1016/j.memsci.2020.118949.[Export]
BibTeX
@article{FUERST2021118949, author = "Fuerst, Thomas F. and Humrickhouse, Paul W. and Taylor, Chase N. and Shimada, Masashi", title = "Surface effects on deuterium permeation through vanadium membranes", journal = "Journal of Membrane Science", volume = "620", pages = "118949", year = "2021", issn = "0376-7388", doi = "10.1016/j.memsci.2020.118949", url = "https://doi.org/10.1016/j.memsci.2020.118949", keywords = "Hydrogen isotopes, Vanadium, Palladium, Fusion technology, Dense metal membranes" }RIS
TY - JOUR AU - Fuerst, Thomas F. AU - Humrickhouse, Paul W. AU - Taylor, Chase N. AU - Shimada, Masashi TI - Surface effects on deuterium permeation through vanadium membranes JO - Journal of Membrane Science PY - 2021 VL - 620 DO - 10.1016/j.memsci.2020.118949 UR - https://doi.org/10.1016/j.memsci.2020.118949 SP - 118949 ER -Plain Text
Thomas F. Fuerst, Paul W. Humrickhouse, Chase N. Taylor, and Masashi Shimada. Surface effects on deuterium permeation through vanadium membranes. Journal of Membrane Science, 620:118949, 2021. URL: https://doi.org/10.1016/j.memsci.2020.118949, doi:10.1016/j.memsci.2020.118949. - S. Fukada, T. Terai, S. Konishi, K. Katayama, T. Chikada, Y. Edao, T. Muroga, M. Shimada, B. Merrill, and D. K. Sze.
Clarification of Tritium Behavior in Pb-Li Blanket System.
Materials Transactions, 54(4):425–429, 2013.
URL: https://doi.org/10.2320/matertrans.MG201203, doi:10.2320/matertrans.MG201203.[Export]
BibTeX
@article{S.Fukada2013MG201203, author = "Fukada, S. and Terai, T. and Konishi, S. and Katayama, K. and Chikada, T. and Edao, Y. and Muroga, T. and Shimada, M. and Merrill, B. and Sze, D. K.", title = "{Clarification of Tritium Behavior in Pb-Li Blanket System}", journal = "Materials Transactions", volume = "54", number = "4", pages = "425--429", year = "2013", doi = "10.2320/matertrans.MG201203", url = "https://doi.org/10.2320/matertrans.MG201203" }RIS
TY - JOUR AU - Fukada, S. AU - Terai, T. AU - Konishi, S. AU - Katayama, K. AU - Chikada, T. AU - Edao, Y. AU - Muroga, T. AU - Shimada, M. AU - Merrill, B. AU - Sze, D. K. TI - Clarification of Tritium Behavior in Pb-Li Blanket System JO - Materials Transactions PY - 2013 VL - 54 IS - 4 DO - 10.2320/matertrans.MG201203 UR - https://doi.org/10.2320/matertrans.MG201203 SP - 425 EP - 429 ER -Plain Text
S. Fukada, T. Terai, S. Konishi, K. Katayama, T. Chikada, Y. Edao, T. Muroga, M. Shimada, B. Merrill, and D. K. Sze. Clarification of Tritium Behavior in Pb-Li Blanket System. Materials Transactions, 54(4):425–429, 2013. URL: https://doi.org/10.2320/matertrans.MG201203, doi:10.2320/matertrans.MG201203. - Paul W. Humrickhouse and Brad J. Merrill.
Vacuum permeator analysis for extraction of tritium from DCLL blankets.
Fusion Science and Technology, 68(2):295–302, 2015.
doi:10.13182/FST14-941.[Export]
BibTeX
@article{Humrickhouse2015PAV, author = "Humrickhouse, Paul W. and Merrill, Brad J.", title = "Vacuum Permeator Analysis for Extraction of Tritium from {DCLL} Blankets", journal = "Fusion Science and Technology", year = "2015", volume = "68", number = "2", pages = "295--302", doi = "10.13182/FST14-941" }RIS
TY - JOUR AU - Humrickhouse, Paul W. AU - Merrill, Brad J. TI - Vacuum Permeator Analysis for Extraction of Tritium from DCLL Blankets JO - Fusion Science and Technology PY - 2015 VL - 68 IS - 2 DO - 10.13182/FST14-941 SP - 295 EP - 302 ER -Plain Text
Paul W. Humrickhouse and Brad J. Merrill. Vacuum permeator analysis for extraction of tritium from DCLL blankets. Fusion Science and Technology, 68(2):295–302, 2015. doi:10.13182/FST14-941. - W. H. Linton, Jr. and T. K. Sherwood.
Mass transfer from solid shapes to water in streamline and turbulent flow.
Chemical Engineering Progress, 46(5):258–264, 1950.[Export]
BibTeX
@article{linton1950mass, author = "Linton, Jr., W. H. and Sherwood, T. K.", title = "Mass Transfer from Solid Shapes to Water in Streamline and Turbulent Flow", journal = "Chemical Engineering Progress", volume = "46", number = "5", pages = "258--264", year = "1950" }RIS
TY - JOUR AU - Linton, W. H., Jr. AU - Sherwood, T. K. TI - Mass Transfer from Solid Shapes to Water in Streamline and Turbulent Flow JO - Chemical Engineering Progress PY - 1950 VL - 46 IS - 5 SP - 258 EP - 264 ER -Plain Text
W. H. Linton, Jr. and T. K. Sherwood. Mass transfer from solid shapes to water in streamline and turbulent flow. Chemical Engineering Progress, 46(5):258–264, 1950. - Ewald Veleckis and Russell Keith Edwards.
Thermodynamic properties in the systems vanadium-hydrogen, niobium-hydrogen, and tantalum-hydrogen.
The Journal of Physical Chemistry, 73(3):683–692, 03 1969.
URL: https://doi.org/10.1021/j100723a033, doi:10.1021/j100723a033.[Export]
BibTeX
@article{10.1021/j100723a033, author = "Veleckis, Ewald and Edwards, Russell Keith", title = "Thermodynamic Properties in the Systems Vanadium-Hydrogen, Niobium-Hydrogen, and Tantalum-Hydrogen", journal = "The Journal of Physical Chemistry", volume = "73", number = "3", pages = "683--692", year = "1969", month = "03", issn = "0022-3654", doi = "10.1021/j100723a033", url = "https://doi.org/10.1021/j100723a033" }RIS
TY - JOUR AU - Veleckis, Ewald AU - Edwards, Russell Keith TI - Thermodynamic Properties in the Systems Vanadium-Hydrogen, Niobium-Hydrogen, and Tantalum-Hydrogen JO - The Journal of Physical Chemistry PY - 1969 VL - 73 IS - 3 DO - 10.1021/j100723a033 UR - https://doi.org/10.1021/j100723a033 SP - 683 EP - 692 ER -Plain Text
Ewald Veleckis and Russell Keith Edwards. Thermodynamic properties in the systems vanadium-hydrogen, niobium-hydrogen, and tantalum-hydrogen. The Journal of Physical Chemistry, 73(3):683–692, 03 1969. URL: https://doi.org/10.1021/j100723a033, doi:10.1021/j100723a033. - J. Völkl and G. Alefeld.
Hydrogen diffusion in metals.
In A. S. Nowick and J. J. Burton, editors, Diffusion in Solids, pages 231–302.
Academic Press, 1975.
URL: https://doi.org/10.1016/B978-0-12-522660-8.50010-6, doi:10.1016/B978-0-12-522660-8.50010-6.[Export]
BibTeX
@incollection{VOLKL1975231, author = {V{\"o}lkl, J. and Alefeld, G.}, editor = "Nowick, A. S. and Burton, J. J.", title = "Hydrogen Diffusion in Metals", booktitle = "Diffusion in Solids", publisher = "Academic Press", pages = "231--302", year = "1975", isbn = "978-0-12-522660-8", doi = "10.1016/B978-0-12-522660-8.50010-6", url = "https://doi.org/10.1016/B978-0-12-522660-8.50010-6" }RIS
TY - CHAP AU - Völkl, J. AU - Alefeld, G. TI - Hydrogen Diffusion in Metals T2 - Diffusion in Solids PY - 1975 PB - Academic Press DO - 10.1016/B978-0-12-522660-8.50010-6 UR - https://doi.org/10.1016/B978-0-12-522660-8.50010-6 SP - 231 EP - 302 ER -Plain Text
J. Völkl and G. Alefeld. Hydrogen diffusion in metals. In A. S. Nowick and J. J. Burton, editors, Diffusion in Solids, pages 231–302. Academic Press, 1975. URL: https://doi.org/10.1016/B978-0-12-522660-8.50010-6, doi:10.1016/B978-0-12-522660-8.50010-6.
(test/tests/ver-1p/ver-1p.i)
# Reference:
# T. F. Fuerst, M. D. Eklund, J. A. Leland, A. A. Riet, C. N. Taylor,
# "Parametric Study of the Vacuum Permeator for the Tritium Extraction
# eXperiment," Fusion Science and Technology 79 (2023) 1224-1234.
# DOI: 10.1080/15361055.2023.2196237
#
# Radial transport:
# bulk PbLi -> PbLi/V partition -> cylindrical V diffusion
# -> recombination -> vacuum
#
# Assumptions:
# - steady state and isothermal operation
# - zero pressure on the vacuum side
# Nominal geometry / hydraulics at 673.15 K and Reynolds number = 1e5:
# Tube length = 5.00 m, inner radius = 4.75 mm,
# outer radius = 5.00 mm, and PbLi velocity = 1.78 m/s.
# ============================================================================
# -----------------------------------------------------------------------------
# Inputs
# -----------------------------------------------------------------------------
temperature = '${units 673.15 K}' # 400 C
gas_constant = '${units 8.31446261815324 J/mol/K}'
concentration_inlet = '${units 1.0 mol/m^3}' # dissolved H-isotope concentration
reynolds_number = 1e5
tube_length = '${units 5.00 m}' # active permeator length
radius_inner_wall = '${units 0.00475 m}'
radius_outer_wall = '${units 0.00500 m}'
number_axial_segments = 20
# -----------------------------------------------------------------------------
# Nominal properties
# -----------------------------------------------------------------------------
diffusivity_PbLi = '${units ${fparse 8.30e-9*exp(-7.37e3/(gas_constant*temperature))} m^2/s}'
density_PbLi = '${units ${fparse 10520.35 - 1.19051*temperature} kg/m^3}'
molar_mass_PbLi_atom = '${units 0.1731558 kg/mol}'
atomic_concentration_PbLi = '${units ${fparse density_PbLi/molar_mass_PbLi_atom} mol/m^3}'
solubility_factor_PbLi = '${units 4.32e-7 1/Pa^0.5}' # Table II in the paper
solubility_PbLi = '${units ${fparse solubility_factor_PbLi*atomic_concentration_PbLi} mol/m^3/Pa^(1/2)}'
diffusivity_vanadium = '${units ${fparse 2.90e-8*exp(-4.20e3/(gas_constant*temperature))} m^2/s}'
solubility_vanadium = '${units ${fparse 0.138*exp(29.0e3/(gas_constant*temperature))} mol/m^3/Pa^(1/2)}'
recombination_coefficient = '${units 3.1582e-9 m^4/mol/s}'
outer_to_inner_area_ratio = ${fparse radius_outer_wall/radius_inner_wall}
recombination_coefficient_inner_area = '${units ${fparse outer_to_inner_area_ratio*recombination_coefficient} m^4/mol/s}'
# -----------------------------------------------------------------------------
# Properties from inputs
# -----------------------------------------------------------------------------
velocity_PbLi = '${units 1.78 m/s}'
hydraulic_diameter = '${units ${fparse 2.0*radius_inner_wall} m}'
flow_area = '${units ${fparse pi*radius_inner_wall^2} m^2}'
volumetric_flow_rate_PbLi = '${units ${fparse velocity_PbLi*flow_area} m^3/s}'
kinematic_viscosity_PbLi = '${units ${fparse velocity_PbLi*hydraulic_diameter/reynolds_number} m^2/s}'
schmidt_number = ${fparse kinematic_viscosity_PbLi/diffusivity_PbLi}
# -----------------------------------------------------------------------------
# PbLi mass transfer from Sherwood correlations
# -----------------------------------------------------------------------------
# sherwood_number = 0.023 reynolds_number^0.83 schmidt_number^(1/3)
# mass_transfer_coefficient_PbLi = sherwood_number diffusivity_PbLi / (2 radius_inner_wall)
sherwood_number = ${fparse 0.023*reynolds_number^0.83*schmidt_number^(1.0/3.0)}
mass_transfer_coefficient_PbLi = '${units ${fparse sherwood_number*diffusivity_PbLi/hydraulic_diameter} m/s}'
# -----------------------------------------------------------------------------
# Radial equations
# -----------------------------------------------------------------------------
#
# Liquid flux on the inner area:
# J_i = mass_transfer_coefficient_PbLi * (C_bulk - C_L2)
#
# Interface partition:
# C_S1 = partition_ratio * C_L2
#
# Cylindrical diffusion expressed on the inner area:
# J_i = diffusivity_vanadium * (C_S1 - C_S2)
# / (radius_inner_wall * ln(radius_outer_wall/radius_inner_wall))
#
# Recombination occurs on the outer vacuum surface:
# J_o = recombination_coefficient * C_S2^2
#
# Conservation of permeation rate requires
# 2*pi*radius_inner_wall*J_i = 2*pi*radius_outer_wall*J_o,
#
# so the equivalent recombination flux on the inner area is
# J_i = recombination_coefficient_inner_area * C_S2^2.
#
# Define
# partition_ratio = solubility_vanadium / solubility_PbLi
# membrane_resistance = radius_inner_wall
# * ln(radius_outer_wall/radius_inner_wall)
# / diffusivity_vanadium
# radial_transport_coefficient =
# partition_ratio / mass_transfer_coefficient_PbLi
# + membrane_resistance
#
# Then the physical positive root is
# C_S2(C) = 2*partition_ratio*C
# / [1 + sqrt(1 + 4*radial_transport_coefficient
# * recombination_coefficient_inner_area
# * partition_ratio*C)].
partition_ratio = ${fparse solubility_vanadium/solubility_PbLi}
membrane_resistance = '${units ${fparse radius_inner_wall*log(radius_outer_wall/radius_inner_wall)/diffusivity_vanadium} s/m}'
radial_transport_coefficient = '${units ${fparse partition_ratio/mass_transfer_coefficient_PbLi+membrane_resistance} s/m}'
# -----------------------------------------------------------------------------
# Control volumes (axial direction)
# -----------------------------------------------------------------------------
axial_segment_length = '${units ${fparse tube_length/number_axial_segments} m}'
segment_permeation_area = '${units ${fparse 2.0*pi*radius_inner_wall*axial_segment_length} m^2}'
total_permeation_area = '${units ${fparse 2.0*pi*radius_inner_wall*tube_length} m^2}'
segment_area_to_flow_ratio = '${units ${fparse segment_permeation_area/volumetric_flow_rate_PbLi} s/m}'
# -----------------------------------------------------------------------------
# Pressure drop
# -----------------------------------------------------------------------------
darcy_friction_factor = ${fparse 0.3164/(reynolds_number^0.25)}
pressure_drop_straight = '${units ${fparse darcy_friction_factor*(tube_length/hydraulic_diameter)*density_PbLi*velocity_PbLi^2/2.0} Pa}'
# -----------------------------------------------------------------------------
# Mesh (Needed just to run the model)
# -----------------------------------------------------------------------------
[Mesh]
type = GeneratedMesh
dim = 1
nx = 1
[]
# -----------------------------------------------------------------------------
# Unknowns (concentrations at the axial direction)
# -----------------------------------------------------------------------------
# concentration_segment_01 is the concentration at z = 0.25 m, ...,
# concentration_segment_20 is the concentration at z = 5.00 m.
[Variables]
[concentration_segment_01]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_02]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_03]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_04]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_05]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_06]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_07]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_08]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_09]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_10]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_11]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_12]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_13]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_14]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_15]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_16]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_17]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_18]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_19]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[concentration_segment_20]
family = SCALAR
initial_condition = ${concentration_inlet}
[]
[]
# -----------------------------------------------------------------------------
# Steady state mass balance
# -----------------------------------------------------------------------------
# For each segment i:
# volumetric_flow_rate_PbLi*C_up = volumetric_flow_rate_PbLi*C_down
# + J(C_up)*segment_permeation_area
# or
# C_down - C_up + segment_area_to_flow_ratio*J(C_up) = 0
#
# The radial flux J(C_up) is evaluated using the local upstream concentration.
[ScalarKernels]
[segment_01]
type = ParsedODEKernel
variable = concentration_segment_01
expression = 'concentration_segment_01 - concentration_upstream + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_upstream/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_upstream)))^2'
constant_names = 'concentration_upstream segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${concentration_inlet} ${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_02]
type = ParsedODEKernel
variable = concentration_segment_02
coupled_variables = 'concentration_segment_01'
expression = 'concentration_segment_02 - concentration_segment_01 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_01/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_01)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_03]
type = ParsedODEKernel
variable = concentration_segment_03
coupled_variables = 'concentration_segment_02'
expression = 'concentration_segment_03 - concentration_segment_02 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_02/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_02)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_04]
type = ParsedODEKernel
variable = concentration_segment_04
coupled_variables = 'concentration_segment_03'
expression = 'concentration_segment_04 - concentration_segment_03 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_03/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_03)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_05]
type = ParsedODEKernel
variable = concentration_segment_05
coupled_variables = 'concentration_segment_04'
expression = 'concentration_segment_05 - concentration_segment_04 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_04/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_04)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_06]
type = ParsedODEKernel
variable = concentration_segment_06
coupled_variables = 'concentration_segment_05'
expression = 'concentration_segment_06 - concentration_segment_05 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_05/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_05)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_07]
type = ParsedODEKernel
variable = concentration_segment_07
coupled_variables = 'concentration_segment_06'
expression = 'concentration_segment_07 - concentration_segment_06 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_06/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_06)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_08]
type = ParsedODEKernel
variable = concentration_segment_08
coupled_variables = 'concentration_segment_07'
expression = 'concentration_segment_08 - concentration_segment_07 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_07/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_07)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_09]
type = ParsedODEKernel
variable = concentration_segment_09
coupled_variables = 'concentration_segment_08'
expression = 'concentration_segment_09 - concentration_segment_08 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_08/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_08)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_10]
type = ParsedODEKernel
variable = concentration_segment_10
coupled_variables = 'concentration_segment_09'
expression = 'concentration_segment_10 - concentration_segment_09 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_09/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_09)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_11]
type = ParsedODEKernel
variable = concentration_segment_11
coupled_variables = 'concentration_segment_10'
expression = 'concentration_segment_11 - concentration_segment_10 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_10/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_10)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_12]
type = ParsedODEKernel
variable = concentration_segment_12
coupled_variables = 'concentration_segment_11'
expression = 'concentration_segment_12 - concentration_segment_11 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_11/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_11)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_13]
type = ParsedODEKernel
variable = concentration_segment_13
coupled_variables = 'concentration_segment_12'
expression = 'concentration_segment_13 - concentration_segment_12 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_12/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_12)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_14]
type = ParsedODEKernel
variable = concentration_segment_14
coupled_variables = 'concentration_segment_13'
expression = 'concentration_segment_14 - concentration_segment_13 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_13/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_13)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_15]
type = ParsedODEKernel
variable = concentration_segment_15
coupled_variables = 'concentration_segment_14'
expression = 'concentration_segment_15 - concentration_segment_14 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_14/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_14)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_16]
type = ParsedODEKernel
variable = concentration_segment_16
coupled_variables = 'concentration_segment_15'
expression = 'concentration_segment_16 - concentration_segment_15 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_15/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_15)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_17]
type = ParsedODEKernel
variable = concentration_segment_17
coupled_variables = 'concentration_segment_16'
expression = 'concentration_segment_17 - concentration_segment_16 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_16/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_16)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_18]
type = ParsedODEKernel
variable = concentration_segment_18
coupled_variables = 'concentration_segment_17'
expression = 'concentration_segment_18 - concentration_segment_17 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_17/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_17)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_19]
type = ParsedODEKernel
variable = concentration_segment_19
coupled_variables = 'concentration_segment_18'
expression = 'concentration_segment_19 - concentration_segment_18 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_18/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_18)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[segment_20]
type = ParsedODEKernel
variable = concentration_segment_20
coupled_variables = 'concentration_segment_19'
expression = 'concentration_segment_20 - concentration_segment_19 + segment_area_to_flow_ratio*recombination_coefficient_inner_area*(2.0*partition_ratio_local*concentration_segment_19/(1.0+sqrt(1.0+4.0*radial_transport_coefficient_local*recombination_coefficient_inner_area*partition_ratio_local*concentration_segment_19)))^2'
constant_names = 'segment_area_to_flow_ratio recombination_coefficient_inner_area partition_ratio_local radial_transport_coefficient_local'
constant_expressions = '${segment_area_to_flow_ratio} ${recombination_coefficient_inner_area} ${partition_ratio} ${radial_transport_coefficient}'
evalerror_behavior = error
[]
[]
# -----------------------------------------------------------------------------
# Postprocessor
# -----------------------------------------------------------------------------
[Postprocessors]
[concentration_inlet_pp]
type = ConstantPostprocessor
value = ${concentration_inlet}
execute_on = 'initial timestep_end'
[]
[concentration_segment_01_pp]
type = ScalarVariable
variable = concentration_segment_01
execute_on = 'initial timestep_end'
[]
[concentration_segment_02_pp]
type = ScalarVariable
variable = concentration_segment_02
execute_on = 'initial timestep_end'
[]
[concentration_segment_03_pp]
type = ScalarVariable
variable = concentration_segment_03
execute_on = 'initial timestep_end'
[]
[concentration_segment_04_pp]
type = ScalarVariable
variable = concentration_segment_04
execute_on = 'initial timestep_end'
[]
[concentration_segment_05_pp]
type = ScalarVariable
variable = concentration_segment_05
execute_on = 'initial timestep_end'
[]
[concentration_segment_06_pp]
type = ScalarVariable
variable = concentration_segment_06
execute_on = 'initial timestep_end'
[]
[concentration_segment_07_pp]
type = ScalarVariable
variable = concentration_segment_07
execute_on = 'initial timestep_end'
[]
[concentration_segment_08_pp]
type = ScalarVariable
variable = concentration_segment_08
execute_on = 'initial timestep_end'
[]
[concentration_segment_09_pp]
type = ScalarVariable
variable = concentration_segment_09
execute_on = 'initial timestep_end'
[]
[concentration_segment_10_pp]
type = ScalarVariable
variable = concentration_segment_10
execute_on = 'initial timestep_end'
[]
[concentration_segment_11_pp]
type = ScalarVariable
variable = concentration_segment_11
execute_on = 'initial timestep_end'
[]
[concentration_segment_12_pp]
type = ScalarVariable
variable = concentration_segment_12
execute_on = 'initial timestep_end'
[]
[concentration_segment_13_pp]
type = ScalarVariable
variable = concentration_segment_13
execute_on = 'initial timestep_end'
[]
[concentration_segment_14_pp]
type = ScalarVariable
variable = concentration_segment_14
execute_on = 'initial timestep_end'
[]
[concentration_segment_15_pp]
type = ScalarVariable
variable = concentration_segment_15
execute_on = 'initial timestep_end'
[]
[concentration_segment_16_pp]
type = ScalarVariable
variable = concentration_segment_16
execute_on = 'initial timestep_end'
[]
[concentration_segment_17_pp]
type = ScalarVariable
variable = concentration_segment_17
execute_on = 'initial timestep_end'
[]
[concentration_segment_18_pp]
type = ScalarVariable
variable = concentration_segment_18
execute_on = 'initial timestep_end'
[]
[concentration_segment_19_pp]
type = ScalarVariable
variable = concentration_segment_19
execute_on = 'initial timestep_end'
[]
[concentration_segment_20_pp]
type = ScalarVariable
variable = concentration_segment_20
execute_on = 'initial timestep_end'
[]
[pressure_drop_straight_Pa]
type = ConstantPostprocessor
value = ${pressure_drop_straight}
execute_on = 'initial timestep_end'
[]
[extraction_efficiency]
type = ParsedPostprocessor
expression = '1.0-concentration_outlet/concentration_inlet_local'
pp_names = 'concentration_segment_20_pp'
pp_symbols = 'concentration_outlet'
constant_names = 'concentration_inlet_local'
constant_expressions = '${concentration_inlet}'
execute_on = 'timestep_end'
[]
[average_permeation_flux_mol_m2_s]
type = ParsedPostprocessor
expression = 'volumetric_flow_rate*(concentration_inlet_local-concentration_outlet)/permeation_area'
pp_names = 'concentration_segment_20_pp'
pp_symbols = 'concentration_outlet'
constant_names = 'volumetric_flow_rate concentration_inlet_local permeation_area'
constant_expressions = '${volumetric_flow_rate_PbLi} ${concentration_inlet} ${total_permeation_area}'
execute_on = 'timestep_end'
[]
[]
[Executioner]
type = Steady
solve_type = NEWTON
nl_abs_tol = 1e-25
nl_rel_tol = 1e-8
nl_max_its = 100
[]
[Outputs]
console = true
[csv]
type = CSV
file_base = ver-1p_out
precision = 16
scientific_notation = true
[]
[]
(test/tests/ver-1p/tests)
[Tests]
design = 'ParsedODEKernel.md ScalarVariable.md ParsedPostprocessor.md'
verification = 'ver-1p.md'
issues = '#457'
[ver-1p_csvdiff]
type = CSVDiff
input = ver-1p.i
csvdiff = 'ver-1p_out.csv'
requirement = 'The system shall solve the steady-state PAV model for PbLi flowing through a vanadium permeator at 673.15 K, including PbLi mass transfer, PbLi/vanadium surface condition, solid diffusion, recombination, and axial concentration depletion.'
[]
[ver-1p_analytical_comparison]
type = RunCommand
command = 'python3 comparison_ver-1p.py'
requirement = 'The system shall compare the 20-segment TMAP8 solution with an independent analytical solution of the PAV radial transport and axial concentration equations.'
required_python_packages = 'matplotlib numpy scipy'
[]
[](test/tests/ver-1p/comparison_ver-1p.py)
#!/usr/bin/env python3
"""Analytical comparison for the Fuerst et al. (2023) PAV test."""
from __future__ import annotations
import csv
import math
import os
import re
import sys
from pathlib import Path
from matplotlib.patches import Rectangle
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import brentq
SCRIPT_FOLDER = Path(os.path.abspath(__file__)).parent
os.chdir(SCRIPT_FOLDER)
_NUMERIC_LITERAL = re.compile(r"^[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?")
def locate_input_file() -> Path:
"""Return the location of the ver-1p input file."""
candidates = [SCRIPT_FOLDER / "ver-1p.i"]
candidates.extend(
parent / "test" / "tests" / "ver-1p" / "ver-1p.i"
for parent in SCRIPT_FOLDER.parents
)
for candidate in candidates:
if candidate.is_file():
return candidate
searched = "\n ".join(str(path) for path in candidates)
raise FileNotFoundError(
f"Could not find ver-1p input file. Searched:\n {searched}"
)
def get_numeric_parameter(parameter_name: str) -> float:
"""Read a numerical top level parameter from ver-1p.i."""
input_file = locate_input_file()
with input_file.open(encoding="utf-8") as handle:
for line in handle:
stripped = line.strip()
if not stripped.startswith(f"{parameter_name} ="):
continue
raw_value = stripped.split("=", maxsplit=1)[1]
raw_value = raw_value.split("#", maxsplit=1)[0].strip()
if (
len(raw_value) >= 2
and raw_value[0] == raw_value[-1]
and raw_value[0] in {"'", '"'}
):
raw_value = raw_value[1:-1].strip()
if raw_value.startswith("${units"):
raw_value = raw_value[len("${units") :].lstrip()
match = _NUMERIC_LITERAL.match(raw_value)
if match is None:
raise ValueError(
f"Parameter {parameter_name!r} in {input_file} is not a "
"direct numerical value that this script can read."
)
return float(match.group(0))
raise KeyError(f"Could not find parameter {parameter_name!r} in {input_file}")
# Model parameters read from ver-1p.i
TEMPERATURE = get_numeric_parameter("temperature")
GAS_CONSTANT = get_numeric_parameter("gas_constant")
CONCENTRATION_INLET = get_numeric_parameter("concentration_inlet")
REYNOLDS_NUMBER = get_numeric_parameter("reynolds_number")
TUBE_LENGTH = get_numeric_parameter("tube_length")
RADIUS_INNER_WALL = get_numeric_parameter("radius_inner_wall")
RADIUS_OUTER_WALL = get_numeric_parameter("radius_outer_wall")
NUMBER_AXIAL_SEGMENTS = int(get_numeric_parameter("number_axial_segments"))
VELOCITY_PBLI = get_numeric_parameter("velocity_PbLi")
SOLUBILITY_FACTOR_PBLI = get_numeric_parameter("solubility_factor_PbLi")
RECOMBINATION_COEFFICIENT = get_numeric_parameter("recombination_coefficient")
MOLAR_MASS_PBLI_ATOM = get_numeric_parameter("molar_mass_PbLi_atom")
OUTER_TO_INNER_AREA_RATIO = RADIUS_OUTER_WALL / RADIUS_INNER_WALL
RECOMBINATION_COEFFICIENT_INNER_AREA = (
OUTER_TO_INNER_AREA_RATIO * RECOMBINATION_COEFFICIENT
)
# Verification criterion for the 20 segment solution
MAX_ERROR_PERCENT = 1.0
def locate_gold_csv() -> Path:
"""Return the location of the ver-1p gold CSV file."""
candidates = [SCRIPT_FOLDER / "gold" / "ver-1p_out.csv"]
candidates.extend(
parent / "test" / "tests" / "ver-1p" / "gold" / "ver-1p_out.csv"
for parent in SCRIPT_FOLDER.parents
)
for candidate in candidates:
if candidate.is_file():
return candidate
searched = "\n ".join(str(path) for path in candidates)
raise FileNotFoundError(f"Could not find ver-1p gold CSV. Searched:\n {searched}")
def surface_concentration(
concentration_bulk: float,
partition_ratio: float,
radial_transport_coefficient: float,
) -> float:
"""Return the vanadium concentration at the vacuum surface."""
radical = (
1.0
+ 4.0
* radial_transport_coefficient
* RECOMBINATION_COEFFICIENT_INNER_AREA
* partition_ratio
* concentration_bulk
)
return 2.0 * partition_ratio * concentration_bulk / (1.0 + math.sqrt(radical))
def analytical_solution() -> dict[str, object]:
"""Return the continuous analytical axial concentration solution."""
# PbLi properties.
diffusivity_pbli = 8.30e-9 * math.exp(-7.37e3 / (GAS_CONSTANT * TEMPERATURE))
density_pbli = 10520.35 - 1.19051 * TEMPERATURE
atomic_concentration_pbli = density_pbli / MOLAR_MASS_PBLI_ATOM
solubility_pbli = SOLUBILITY_FACTOR_PBLI * atomic_concentration_pbli
# Vanadium properties.
diffusivity_vanadium = 2.90e-8 * math.exp(-4.2e3 / (GAS_CONSTANT * TEMPERATURE))
solubility_vanadium = 0.138 * math.exp(29.0e3 / (GAS_CONSTANT * TEMPERATURE))
# Liquid mass transfer.
hydraulic_diameter = 2.0 * RADIUS_INNER_WALL
volumetric_flow_rate_pbli = VELOCITY_PBLI * math.pi * RADIUS_INNER_WALL**2
kinematic_viscosity_pbli = VELOCITY_PBLI * hydraulic_diameter / REYNOLDS_NUMBER
schmidt_number = kinematic_viscosity_pbli / diffusivity_pbli
sherwood_number = 0.023 * REYNOLDS_NUMBER**0.83 * schmidt_number ** (1.0 / 3.0)
mass_transfer_coefficient_pbli = (
sherwood_number * diffusivity_pbli / hydraulic_diameter
)
# Radial transport parameters.
partition_ratio = solubility_vanadium / solubility_pbli
membrane_resistance = (
RADIUS_INNER_WALL
* math.log(RADIUS_OUTER_WALL / RADIUS_INNER_WALL)
/ diffusivity_vanadium
)
radial_transport_coefficient = (
partition_ratio / mass_transfer_coefficient_pbli + membrane_resistance
)
surface_concentration_inlet = surface_concentration(
CONCENTRATION_INLET,
partition_ratio,
radial_transport_coefficient,
)
def axial_residual(
surface_concentration_value: float,
axial_position: float,
) -> float:
return (
1.0 / surface_concentration_value
- 1.0 / surface_concentration_inlet
- 2.0
* radial_transport_coefficient
* RECOMBINATION_COEFFICIENT_INNER_AREA
* math.log(surface_concentration_value / surface_concentration_inlet)
- (
2.0
* math.pi
* RADIUS_INNER_WALL
* partition_ratio
* RECOMBINATION_COEFFICIENT_INNER_AREA
/ volumetric_flow_rate_pbli
)
* axial_position
)
concentrations = []
for index in range(1, NUMBER_AXIAL_SEGMENTS + 1):
axial_position = index * TUBE_LENGTH / NUMBER_AXIAL_SEGMENTS
surface_concentration_value = brentq(
axial_residual,
surface_concentration_inlet * 1.0e-8,
surface_concentration_inlet,
args=(axial_position,),
xtol=1.0e-14,
rtol=1.0e-14,
)
concentration_bulk = (
surface_concentration_value
+ radial_transport_coefficient
* RECOMBINATION_COEFFICIENT_INNER_AREA
* surface_concentration_value**2
) / partition_ratio
concentrations.append(concentration_bulk)
concentration_outlet = concentrations[-1]
efficiency = 1.0 - concentration_outlet / CONCENTRATION_INLET
# Average flux on the inner surface.
total_permeation_area = 2.0 * math.pi * RADIUS_INNER_WALL * TUBE_LENGTH
average_permeation_flux = (
volumetric_flow_rate_pbli
* (CONCENTRATION_INLET - concentration_outlet)
/ total_permeation_area
)
return {
"concentrations": concentrations,
"outlet": concentration_outlet,
"efficiency": efficiency,
"average_flux": average_permeation_flux,
}
def read_last_csv_row(path: Path) -> dict[str, str]:
"""Return the final data row in a CSV file."""
with path.open(newline="", encoding="utf-8") as stream:
rows = list(csv.DictReader(stream))
if not rows:
raise RuntimeError(f"No data rows found in {path}")
return rows[-1]
def value(row: dict[str, str], name: str) -> float:
"""Return a named CSV value as a float."""
if name not in row:
raise KeyError(f"Column {name!r} is missing from the TMAP8 CSV")
return float(row[name])
def relative_error_percent(
actual: float,
expected: float,
) -> float:
"""Return the absolute relative error as a percentage."""
return abs((actual - expected) / expected) * 100.0
def make_schematic() -> None:
"""Generate a schematic of the segmented PAV verification model."""
figure = plt.figure(figsize=(11.0, 6.4))
axis = figure.add_axes((0.03, 0.04, 0.94, 0.92))
axis.set_xlim(0.0, 12.0)
axis.set_ylim(0.0, 7.0)
axis.axis("off")
highlight_color = "#dce8f5"
highlight_edge = "#6f8fae"
membrane_color = "#d9d9d9"
vacuum_color = "#eef5e9"
radial_arrow_color = "#b94a48"
# ------------------------------------------------------------------
# Vertical axial permeator
# ------------------------------------------------------------------
pipe_left = 2.55
pipe_bottom = 0.85
pipe_width = 1.00
pipe_height = 5.20
pipe_right = pipe_left + pipe_width
pipe_top = pipe_bottom + pipe_height
pipe_center = pipe_left + pipe_width / 2.0
segment_height = pipe_height / NUMBER_AXIAL_SEGMENTS
axis.add_patch(
Rectangle(
(pipe_left, pipe_bottom),
pipe_width,
pipe_height,
facecolor="white",
edgecolor="black",
linewidth=1.8,
zorder=1,
)
)
for segment_index in range(1, NUMBER_AXIAL_SEGMENTS):
y_position = pipe_bottom + segment_index * segment_height
axis.plot(
[pipe_left, pipe_right],
[y_position, y_position],
color="0.70",
linewidth=0.6,
zorder=2,
)
# ------------------------------------------------------------------
# Axial PbLi flow
# ------------------------------------------------------------------
axis.annotate(
"",
xy=(pipe_center, pipe_top - 0.16),
xytext=(pipe_center, pipe_bottom + 0.16),
arrowprops={
"arrowstyle": "->",
"linewidth": 2.0,
"color": "black",
},
zorder=4,
)
axis.text(
pipe_center,
pipe_bottom - 0.48,
"PbLi + T",
ha="center",
va="top",
fontsize=11,
)
axis.text(
pipe_left - 0.15,
pipe_bottom,
r"$z=0$",
ha="right",
va="center",
fontsize=10,
)
axis.text(
pipe_left - 0.15,
pipe_top,
rf"$z={TUBE_LENGTH:.0f}\ \mathrm{{m}}$",
ha="right",
va="center",
fontsize=10,
)
axis.text(
pipe_right + 0.16,
pipe_bottom,
r"$C_{\mathrm{in}}$",
ha="left",
va="center",
fontsize=12,
)
axis.text(
pipe_right + 0.16,
pipe_top,
r"$C_{\mathrm{out}}$",
ha="left",
va="center",
fontsize=12,
)
axis.text(
pipe_center,
pipe_top + 0.45,
"20 axial control volumes",
ha="center",
fontsize=12,
fontweight="bold",
)
# ------------------------------------------------------------------
# Show one representative axial control volume
# ------------------------------------------------------------------
selected_segment = NUMBER_AXIAL_SEGMENTS // 2
selected_bottom = pipe_bottom + selected_segment * segment_height
selected_top = selected_bottom + segment_height
selected_center = (selected_bottom + selected_top) / 2.0
highlight_padding_x = 0.12
highlight_padding_y = 0.07
highlighted_left = pipe_left - highlight_padding_x
highlighted_right = pipe_right + highlight_padding_x
highlighted_bottom = selected_bottom - highlight_padding_y
highlighted_top = selected_top + highlight_padding_y
axis.add_patch(
Rectangle(
(highlighted_left, highlighted_bottom),
highlighted_right - highlighted_left,
highlighted_top - highlighted_bottom,
facecolor=highlight_color,
edgecolor=highlight_edge,
linewidth=2.2,
zorder=3,
)
)
# Redraw axial flow through the highlighted segment.
axis.annotate(
"",
xy=(pipe_center, selected_top - 0.02),
xytext=(pipe_center, selected_bottom + 0.02),
arrowprops={
"arrowstyle": "->",
"linewidth": 1.5,
"color": "black",
},
zorder=5,
)
# ------------------------------------------------------------------
# Axial segment balance
# ------------------------------------------------------------------
relationship_x = 1.18
axis.text(
relationship_x,
selected_center + 0.68,
"Axial segment balance",
ha="center",
va="center",
fontsize=10,
fontweight="bold",
)
axis.text(
relationship_x,
selected_center + 0.20,
(r"$C_i-C_{i-1}" r"+\frac{A_{\mathrm{seg}}}{Q}" r"J(C_{i-1})=0$"),
ha="center",
va="center",
fontsize=10,
)
axis.text(
relationship_x,
selected_center - 0.45,
"Outlet of segment $i-1$\n" "is the inlet of segment $i$",
ha="center",
va="center",
fontsize=8.5,
)
segment_length = TUBE_LENGTH / NUMBER_AXIAL_SEGMENTS
# ------------------------------------------------------------------
# Symmetric radial permeation from the selected axial segment
# ------------------------------------------------------------------
radial_arrow_y = (
selected_center - 0.07,
selected_center,
selected_center + 0.07,
)
for y_position in radial_arrow_y:
axis.annotate(
"",
xy=(highlighted_left - 0.18, y_position),
xytext=(highlighted_left + 0.05, y_position),
arrowprops={
"arrowstyle": "->",
"linewidth": 1.3,
"color": radial_arrow_color,
},
zorder=5,
)
axis.annotate(
"",
xy=(highlighted_right + 0.18, y_position),
xytext=(highlighted_right - 0.05, y_position),
arrowprops={
"arrowstyle": "->",
"linewidth": 1.3,
"color": radial_arrow_color,
},
zorder=5,
)
model_left = 4.70
model_right = 11.60
model_bottom = 0.85
model_top = 5.75
axis.add_patch(
Rectangle(
(model_left, model_bottom),
model_right - model_left,
model_top - model_bottom,
facecolor=highlight_color,
edgecolor=highlight_edge,
linewidth=1.8,
zorder=-3,
)
)
# Connector lines from the highlighted axial segment.
axis.plot(
[highlighted_right, model_left],
[highlighted_top, model_top - 0.40],
color=highlight_edge,
linewidth=1.2,
)
axis.plot(
[highlighted_right, model_left],
[highlighted_bottom, model_bottom + 0.40],
color=highlight_edge,
linewidth=1.2,
)
axis.text(
(model_left + model_right) / 2.0,
model_top - 0.28,
"Radial transport within one axial control volume",
ha="center",
va="center",
fontsize=11,
fontweight="bold",
)
# ------------------------------------------------------------------
# Radial regions
# ------------------------------------------------------------------
pbli_left = 5.05
radius_inner_location = 7.75
radius_outer_location = 9.65
vacuum_right = 11.25
region_bottom = 1.20
region_top = 4.95
axis.add_patch(
Rectangle(
(pbli_left, region_bottom),
radius_inner_location - pbli_left,
region_top - region_bottom,
facecolor=highlight_color,
edgecolor="none",
zorder=-2,
)
)
axis.add_patch(
Rectangle(
(radius_inner_location, region_bottom),
radius_outer_location - radius_inner_location,
region_top - region_bottom,
facecolor=membrane_color,
edgecolor="none",
zorder=-2,
)
)
axis.add_patch(
Rectangle(
(radius_outer_location, region_bottom),
vacuum_right - radius_outer_location,
region_top - region_bottom,
facecolor=vacuum_color,
edgecolor="none",
zorder=-2,
)
)
# Interface boundaries.
axis.plot(
[radius_inner_location, radius_inner_location],
[region_bottom, region_top],
color="black",
linewidth=1.3,
)
axis.plot(
[radius_outer_location, radius_outer_location],
[region_bottom, region_top],
color="black",
linewidth=1.3,
)
# Region centers.
pbli_center = (pbli_left + radius_inner_location) / 2.0
membrane_center = (radius_inner_location + radius_outer_location) / 2.0
vacuum_center = (radius_outer_location + vacuum_right) / 2.0
# ------------------------------------------------------------------
# Region titles
# ------------------------------------------------------------------
axis.text(
pbli_center,
4.62,
"Bulk PbLi",
ha="center",
fontsize=11,
fontweight="bold",
)
axis.text(
membrane_center,
4.62,
"Membrane",
ha="center",
fontsize=11,
fontweight="bold",
)
axis.text(
vacuum_center,
4.62,
"Vacuum",
ha="center",
fontsize=11,
fontweight="bold",
)
# ------------------------------------------------------------------
# Process descriptions
# ------------------------------------------------------------------
axis.text(
pbli_center,
3.88,
"Liquid mass transfer",
ha="center",
fontsize=9,
)
axis.text(
membrane_center,
3.88,
"Cylindrical diffusion",
ha="center",
fontsize=9,
)
transport_y = 3.22
axis.annotate(
"",
xy=(10.95, transport_y),
xytext=(5.30, transport_y),
arrowprops={
"arrowstyle": "->",
"linewidth": 2.0,
"color": "black",
},
)
# ------------------------------------------------------------------
# Transport equations
# ------------------------------------------------------------------
axis.text(
pbli_center,
2.62,
r"$J_i=K_T(C-C_{L2})$",
ha="center",
fontsize=10,
)
axis.text(
membrane_center,
2.62,
r"$C_S$",
ha="center",
fontsize=11,
)
axis.text(
vacuum_center,
2.62,
r"$J_o=K_R C_{S2}^{\,2}$",
ha="center",
fontsize=10,
)
# ------------------------------------------------------------------
# Interface information
axis.text(
radius_inner_location - 0.28,
2.04,
"PbLi / V interface",
ha="right",
va="center",
fontsize=9,
)
axis.text(
radius_inner_location - 0.28,
1.66,
r"$C_{S1}=(K_S/K_L)C_{L2}$",
ha="right",
va="center",
fontsize=9,
)
axis.text(
radius_inner_location - 0.18,
1.32,
r"$r=r_i$",
ha="right",
fontsize=10,
)
axis.text(
radius_outer_location + 0.28,
2.04,
"Vacuum surface",
ha="left",
va="center",
fontsize=9,
)
axis.text(
vacuum_center,
1.66,
r"$p_{\mathrm{vac}}=0$",
ha="center",
va="center",
fontsize=9,
)
axis.text(
radius_outer_location + 0.18,
1.32,
r"$r=r_o$",
ha="left",
fontsize=10,
)
figure.savefig(
"ver-1p_schematic.png",
dpi=300,
bbox_inches="tight",
)
plt.close(figure)
def make_plot(
actual: dict[str, str],
expected: dict[str, object],
) -> float:
"""Plot TMAP8 and analytical concentration profiles."""
axial_positions = np.linspace(
0.0,
TUBE_LENGTH,
NUMBER_AXIAL_SEGMENTS + 1,
)
tmap8_concentrations = np.array(
[value(actual, "concentration_inlet_pp")]
+ [
value(actual, f"concentration_segment_{index:02d}_pp")
for index in range(1, NUMBER_AXIAL_SEGMENTS + 1)
]
)
analytical_concentrations = np.array(
[CONCENTRATION_INLET] + list(expected["concentrations"])
)
relative_difference = (
tmap8_concentrations - analytical_concentrations
) / analytical_concentrations
rmspe = float(np.sqrt(np.mean(relative_difference**2)) * 100.0)
plt.figure(figsize=(8, 5.4))
plt.plot(
axial_positions,
tmap8_concentrations,
label=f"TMAP8 ({NUMBER_AXIAL_SEGMENTS} segments)",
linewidth=2.2,
)
plt.plot(
axial_positions,
analytical_concentrations,
"--",
label="Analytical",
linewidth=2.0,
)
plt.xlabel("Axial location (m)")
plt.xlim(0.0, TUBE_LENGTH)
plt.ylabel(r"Bulk concentration (mol m$^{-3}$)")
plt.grid(
which="major",
linestyle="--",
alpha=0.3,
)
plt.minorticks_on()
plt.legend()
plt.text(
0.05,
0.08,
f"RMSPE = {rmspe:.3f}%",
transform=plt.gca().transAxes,
)
plt.tight_layout()
plt.savefig(
"ver-1p_comparison_analytical_concentration.png",
dpi=300,
)
plt.close()
return rmspe
def main() -> int:
"""Run the analytical comparison and return an exit status."""
try:
gold_csv = locate_gold_csv()
analytical = analytical_solution()
actual = read_last_csv_row(gold_csv)
except (
FileNotFoundError,
KeyError,
OSError,
RuntimeError,
ValueError,
) as error:
print(f"ERROR: {error}", file=sys.stderr)
return 2
rmspe = make_plot(actual, analytical)
make_schematic()
tmap8_outlet = value(actual, "concentration_segment_20_pp")
tmap8_efficiency = value(actual, "extraction_efficiency")
tmap8_flux = value(actual, "average_permeation_flux_mol_m2_s")
outlet_error = relative_error_percent(
tmap8_outlet,
analytical["outlet"],
)
efficiency_error = relative_error_percent(
tmap8_efficiency,
analytical["efficiency"],
)
flux_error = relative_error_percent(
tmap8_flux,
analytical["average_flux"],
)
print(f"Gold CSV: {gold_csv}")
print("Analytical outlet concentration: " f"{analytical['outlet']:.12e} mol/m^3")
print(f"TMAP8 outlet concentration: {tmap8_outlet:.12e} mol/m^3")
print(f"Outlet relative error: {outlet_error:.6f}%")
print("Analytical extraction efficiency: " f"{analytical['efficiency']:.12e}")
print(f"TMAP8 extraction efficiency: {tmap8_efficiency:.12e}")
print(f"Efficiency relative error: {efficiency_error:.6f}%")
print(
"Analytical average permeation flux: "
f"{analytical['average_flux']:.12e} mol/m^2/s"
)
print("TMAP8 average permeation flux: " f"{tmap8_flux:.12e} mol/m^2/s")
print(f"Flux relative error: {flux_error:.6f}%")
print(f"Concentration profile RMSPE: {rmspe:.6f}%")
failures = []
if rmspe > MAX_ERROR_PERCENT:
failures.append(f"profile RMSPE = {rmspe:.6f}%")
if outlet_error > MAX_ERROR_PERCENT:
failures.append(f"outlet error = {outlet_error:.6f}%")
if efficiency_error > MAX_ERROR_PERCENT:
failures.append(f"efficiency error = {efficiency_error:.6f}%")
if flux_error > MAX_ERROR_PERCENT:
failures.append(f"average flux error = {flux_error:.6f}%")
if failures:
print("\nFAILED analytical comparison:", file=sys.stderr)
for failure in failures:
print(f" {failure}", file=sys.stderr)
print(
f"Allowed error = {MAX_ERROR_PERCENT:.2f}%",
file=sys.stderr,
)
return 1
print(
"TMAP8 agrees with the continuous analytical "
f"solution within {MAX_ERROR_PERCENT:.1f}%."
)
return 0
if __name__ == "__main__":
sys.exit(main())