From b1010082041be77702f8a059cfeb257b0a659f83 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 27 Oct 2025 17:16:41 -0400 Subject: [PATCH 01/26] added GLC model but needs some verification --- src/python/custom_pathsim_blocks.py | 232 ++++++++++++++++++++++++++++ src/python/pathsim_utils.py | 10 +- 2 files changed, 238 insertions(+), 4 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 17820e3..71ddcf8 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -402,3 +402,235 @@ def update(self, t): self.outputs["flux_0"] = flux_0 self.outputs["flux_L"] = flux_L return self.outputs + + +# GLC block + +import numpy as np +from scipy.integrate import solve_bvp + +R = 8.314 # J/(mol·K), Universal gas constant + + +def solve(params): + def solve_tritium_extraction(dimensionless_params, y_T2_in): + """ + Solves the BVP for tritium extraction in a bubble column. + + Args: + params (dict): A dictionary containing the dimensionless parameters: + Bo_l, phi_l, Bo_g, phi_g, psi, nu. + y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-). + + Returns: + sol: The solution object from scipy.integrate.solve_bvp. + """ + + Bo_l = dimensionless_params["Bo_l"] + phi_l = dimensionless_params["phi_l"] + Bo_g = dimensionless_params["Bo_g"] + phi_g = dimensionless_params["phi_g"] + psi = dimensionless_params["psi"] + nu = dimensionless_params["nu"] + + def ode_system(xi, S): + """ + Defines the system of 4 first-order ODEs. + S[0] = x_T (dimensionless liquid concentration) + S[1] = dx_T/d(xi) + S[2] = y_T2 (dimensionless gas concentration) + S[3] = dy_T2/d(xi) + """ + x_T, dx_T_dxi, y_T2, dy_T2_dxi = S + + # Dimensionless driving force theta. This is equation (12) in the paper. + theta = x_T - np.sqrt(np.maximum(0, (1 - psi * xi) * y_T2 / nu)) + + # Equation for d(S[0])/d(xi) = d(x_T)/d(xi) + dS0_dxi = dx_T_dxi + + # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2 from eq (10) + dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) + + # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi) + dS2_dxi = dy_T2_dxi + + # Equation for d(S[3])/d(xi) = d^2(y_T2)/d(xi)^2 from eq (11) + # Avoid division by zero if (1 - psi * xi) is close to zero at xi=1 + denominator = 1 - psi * xi + denominator = np.where(np.isclose(denominator, 0), 1e-9, denominator) + + term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi + term2 = phi_g * theta + dS3_dxi = (Bo_g / denominator) * (term1 - term2) + + return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi)) + + def boundary_conditions(Sa, Sb): + """ + Defines the boundary conditions for the BVP. + Sa: solution at xi = 0 (liquid outlet) + Sb: solution at xi = 1 (liquid inlet) + """ + # Residuals that should be zero for a valid solution. + # Based on equations (16) and (17) in the paper. + + # At xi = 0: dx_T/d(xi) = 0 + res1 = Sa[1] + + # At xi = 1: x_T(1) = 1 - (1/Bo_l) * dx_T/d(xi)|_1 + res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1]) + + # At xi = 0: y_T2(0) = y_T2(0-) + (1/Bo_g) * dy_T2/d(xi)|_0 + res3 = Sa[2] - (y_T2_in + (1 / Bo_g) * Sa[3]) + + # At xi = 1: dy_T2/d(xi) = 0 + res4 = Sb[3] + + return np.array([res1, res2, res3, res4]) + + # Set up the mesh and an initial guess for the solver. + xi = np.linspace(0, 1, 51) + + # A flat initial guess is often good enough to get the solver started. + y_guess = np.zeros((4, xi.size)) + y_guess[0, :] = 0.8 # Guess for x_T + y_guess[2, :] = 1e-5 # Guess for y_T2 + + # Run the BVP solver + sol = solve_bvp(ode_system, boundary_conditions, xi, y_guess, tol=1e-5) + + return sol + + # Unpack parameters + c_T_inlet = params["c_T_inlet"] + y_T2_in = params["y_T2_in"] + P_0 = params["P_0"] + L = params["L"] + u_l = params["u_l"] + u_g0 = params["u_g0"] + ε_g = params["ε_g"] + ε_l = params["ε_l"] + E_g = params["E_g"] + E_l = params["E_l"] + a = params["a"] + h_l = params["h_l"] + ρ_l = params["ρ_l"] + K_s = params["K_s"] + g = params["g"] + D = params["D"] + T = params["T"] + + # Calculate some other parameters + + A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column + A_l = A * ε_l # m^2, Cross-sectional area of the liquid phase + A_g = A * ε_g # m^2, Cross-sectional area of the gas phase + + Q_l = u_l * A_l # m^3/s, Volumetric flow rate of liquid phase + + # Calculate dimensionless values + # Eq 13 in the paper + Bo_l = u_l * L / ((1 - ε_g) * E_l) # Bodenstein number, liquid phase + Φ_l = a * h_l * L / u_l # Mass transfer parameter, liquid phase + # Eq 14 in the paper + Bo_g = ( + u_g0 * L / (ε_g * E_g) + ) # Bodenstein number, gas phase (assumed near plug-flow condition) + Φ_g = ( + 1 / 2 * (R * T * c_T_inlet / P_0) * a * h_l * L / u_g0 + ) # Mass transfer parameter, gas phase + # Eq 15 in the paper + ψ = (ρ_l * g * (1 - ε_g) * L) / P_0 # Hydrostatic pressure ratio + ν = (c_T_inlet / K_s) ** 2 / P_0 # Tritium partial pressure ratio + + dimensionless_params = { + "Bo_l": Bo_l, + "phi_l": Φ_l, + "Bo_g": Bo_g, + "phi_g": Φ_g, + "psi": ψ, + "nu": ν, + } + + # Solve the model + solution = solve_tritium_extraction(dimensionless_params, y_T2_in) + + # --- Results --- + if solution.success: + # --- Dimensionless Results --- + x_T_outlet_dimless = solution.y[0, 0] + efficiency = 1 - x_T_outlet_dimless + y_T2_outlet_gas = solution.y[2, -1] # y_T2 at xi=1 + + # --- Dimensional Results --- + # Liquid concentration at outlet (xi=0) + c_T_outlet = x_T_outlet_dimless * c_T_inlet + + # Gas partial pressure at outlet (xi=1) + P_outlet_gas = P_0 * (1 - dimensionless_params["psi"]) + P_T2_outlet_gas = y_T2_outlet_gas * P_outlet_gas + + # Tritium extraction rate (mol/s) + extraction_rate = Q_l * (c_T_inlet - c_T_outlet) # mol/s + + # Gas velocity at outlet (xi=1) + u_g_outlet = u_g0 / (1 - dimensionless_params["psi"]) + + results = { + "extraction_efficiency [%]": efficiency * 100, + "extraction_rate [mol/s]": extraction_rate, + "c_T_outlet [mol/m^3]": c_T_outlet, + "liquid_vol_flow [m^3/s]": Q_l, + "total_gas_P_outlet [Pa]": P_outlet_gas, + "P_T2_outlet_gas [Pa]": P_T2_outlet_gas, + "gas_vol_flow_outlet [m^3/s]": u_g_outlet * A_g, + } + + else: + print("BVP solver failed to converge.") + + return results + + +from pathsim.blocks import Function + + +class GLC(Function): + """ + This is the docs + """ + + def __init__( + self, P_0, L, u_l, u_g0, ε_g, ε_l, E_g, E_l, a, h_l, ρ_l, K_s, D, T, g=9.81 + ): + self.params = { + "P_0": P_0, + "L": L, + "u_l": u_l, + "u_g0": u_g0, + "ε_g": ε_g, + "ε_l": ε_l, + "E_g": E_g, + "E_l": E_l, + "a": a, + "h_l": h_l, + "ρ_l": ρ_l, + "K_s": K_s, + "g": g, + "D": D, + "T": T, + } + super().__init__(func=self.func) + + def func(self, c_T_inlet, y_T2_inlet): + new_params = self.params.copy() + new_params["c_T_inlet"] = c_T_inlet + new_params["y_T2_in"] = y_T2_inlet + + res = solve(new_params) + + c_T_outlet = res["c_T_outlet [mol/m^3]"] + P_T2_outlet = res["P_T2_outlet_gas [Pa]"] + + return c_T_outlet, P_T2_outlet diff --git a/src/python/pathsim_utils.py b/src/python/pathsim_utils.py index c8ed281..9e9525c 100644 --- a/src/python/pathsim_utils.py +++ b/src/python/pathsim_utils.py @@ -53,6 +53,7 @@ Bubbler, FestimWall, Integrator, + GLC, ) import inspect @@ -134,6 +135,7 @@ "fir": pathsim.blocks.FIR, "interface": pathsim.subsystem.Interface, "switch": pathsim.blocks.Switch, + "glc": GLC, } math_blocks = { @@ -517,12 +519,12 @@ def get_input_index(block: Block, edge: dict, block_to_input_index: dict) -> int return edge["targetHandle"] # TODO maybe we could directly use the targetHandle as a port alias for these: - if type(block) in (Function, ODE, pathsim.blocks.Switch): + if type(block) in (Function, ODE, pathsim.blocks.Switch, GLC): return int(edge["targetHandle"].replace("target-", "")) if isinstance(block, Adder): if block.operations: return int(edge["targetHandle"].replace("target-", "")) - + # make sure that the target block has only one input port (ie. that targetHandle is None) assert edge["targetHandle"] is None, ( f"Target block {block.id} has multiple input ports, " @@ -562,11 +564,11 @@ def get_output_index(block: Block, edge: dict) -> int: ) return output_index # TODO maybe we could directly use the targetHandle as a port alias for these: - if type(block) in (Function, ODE): + if type(block) in (Function, ODE, GLC): # Function and ODE outputs are always in order, so we can use the handle directly assert edge["sourceHandle"], edge return int(edge["sourceHandle"].replace("source-", "")) - + # make sure that the source block has only one output port (ie. that sourceHandle is None) assert edge["sourceHandle"] is None, ( f"Source block {block.id} has multiple output ports, " From 528230642f0b023527b2ed70922b47efa550e5fd Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Oct 2025 13:21:07 -0400 Subject: [PATCH 02/26] added glc to nodeconfig --- src/nodeConfig.js | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/nodeConfig.js b/src/nodeConfig.js index d2e8a8f..4dc91d9 100644 --- a/src/nodeConfig.js +++ b/src/nodeConfig.js @@ -64,6 +64,7 @@ export const nodeTypes = { ode: DynamicHandleNode, interface: DynamicHandleNode, switch: SwitchNode, + glc: createFunctionNode(2, 2), }; export const nodeMathTypes = { @@ -117,7 +118,7 @@ export const nodeCategories = { description: 'Filter and flow control nodes' }, 'Fuel Cycle': { - nodes: ['process', 'process_horizontal', 'bubbler', 'wall'], + nodes: ['process', 'process_horizontal', 'bubbler', 'wall', 'glc'], description: 'Fuel cycle specific nodes' }, 'Others': { @@ -194,6 +195,7 @@ export const getNodeDisplayName = (nodeType) => { 'samplehold': 'Sample Hold', 'comparator': 'Comparator', 'interface': 'Interface', + 'glc': 'Gas-Liquid Contactor', }; return displayNames[nodeType] || nodeType.charAt(0).toUpperCase() + nodeType.slice(1); From 5afac25f50946c175eac58cc89ce3322b5e82d78 Mon Sep 17 00:00:00 2001 From: rossmacdonald98 Date: Tue, 28 Oct 2025 17:34:48 +0000 Subject: [PATCH 03/26] fixes for glc --- src/python/custom_pathsim_blocks.py | 135 ++++++++++++++++++---------- 1 file changed, 86 insertions(+), 49 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 71ddcf8..cb7ff12 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -408,11 +408,13 @@ def update(self, t): import numpy as np from scipy.integrate import solve_bvp +from scipy import constants as const R = 8.314 # J/(mol·K), Universal gas constant def solve(params): + def solve_tritium_extraction(dimensionless_params, y_T2_in): """ Solves the BVP for tritium extraction in a bubble column. @@ -443,13 +445,13 @@ def ode_system(xi, S): """ x_T, dx_T_dxi, y_T2, dy_T2_dxi = S - # Dimensionless driving force theta. This is equation (12) in the paper. + # Dimensionless driving force theta. Eq. 8.8 theta = x_T - np.sqrt(np.maximum(0, (1 - psi * xi) * y_T2 / nu)) # Equation for d(S[0])/d(xi) = d(x_T)/d(xi) dS0_dxi = dx_T_dxi - # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2 from eq (10) + # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2 dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi) @@ -460,9 +462,9 @@ def ode_system(xi, S): denominator = 1 - psi * xi denominator = np.where(np.isclose(denominator, 0), 1e-9, denominator) - term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi - term2 = phi_g * theta - dS3_dxi = (Bo_g / denominator) * (term1 - term2) + term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi # Part of Eq. 9.3.3 (fourth line) + term2 = phi_g * theta # Part of Eq. 9.3.3 (fourth line) + dS3_dxi = (Bo_g / denominator) * (term1 - term2) # Eq. 9.3.3 (fourth line) return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi)) @@ -476,16 +478,16 @@ def boundary_conditions(Sa, Sb): # Based on equations (16) and (17) in the paper. # At xi = 0: dx_T/d(xi) = 0 - res1 = Sa[1] + res1 = Sa[1] # Eq. 10.1 # At xi = 1: x_T(1) = 1 - (1/Bo_l) * dx_T/d(xi)|_1 - res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1]) + res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1]) # Eq. 10.2 # At xi = 0: y_T2(0) = y_T2(0-) + (1/Bo_g) * dy_T2/d(xi)|_0 - res3 = Sa[2] - (y_T2_in + (1 / Bo_g) * Sa[3]) + res3 = Sa[2] - (y_T2_in + (1 / Bo_g) * Sa[3]) # Eq. 10.3 # At xi = 1: dy_T2/d(xi) = 0 - res4 = Sb[3] + res4 = Sb[3] # Eq. 10.4 return np.array([res1, res2, res3, res4]) @@ -504,51 +506,53 @@ def boundary_conditions(Sa, Sb): # Unpack parameters c_T_inlet = params["c_T_inlet"] - y_T2_in = params["y_T2_in"] + P_T2_in = params["P_T2_in"] P_0 = params["P_0"] + ρ_l = params["ρ_l"] + K_s = params["K_s"] + L = params["L"] - u_l = params["u_l"] - u_g0 = params["u_g0"] + D = params["D"] ε_g = params["ε_g"] - ε_l = params["ε_l"] + + Q_l = params["Q_l"] + Q_g = params["Q_g"] + + a = params["a"] + E_g = params["E_g"] E_l = params["E_l"] - a = params["a"] + h_l = params["h_l"] - ρ_l = params["ρ_l"] - K_s = params["K_s"] + g = params["g"] - D = params["D"] T = params["T"] - # Calculate some other parameters + R = params.get("R", R) # Ideal gas constant, J/(mol.K) + # Calculate the superficial flow velocities A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column - A_l = A * ε_l # m^2, Cross-sectional area of the liquid phase - A_g = A * ε_g # m^2, Cross-sectional area of the gas phase - - Q_l = u_l * A_l # m^3/s, Volumetric flow rate of liquid phase + u_g0 = Q_g / A # m/s, superficial gas inlet velocity + u_l = Q_l / A # m/s, superficial liquid inlet velocity # Calculate dimensionless values - # Eq 13 in the paper - Bo_l = u_l * L / ((1 - ε_g) * E_l) # Bodenstein number, liquid phase - Φ_l = a * h_l * L / u_l # Mass transfer parameter, liquid phase - # Eq 14 in the paper - Bo_g = ( - u_g0 * L / (ε_g * E_g) - ) # Bodenstein number, gas phase (assumed near plug-flow condition) - Φ_g = ( - 1 / 2 * (R * T * c_T_inlet / P_0) * a * h_l * L / u_g0 - ) # Mass transfer parameter, gas phase - # Eq 15 in the paper - ψ = (ρ_l * g * (1 - ε_g) * L) / P_0 # Hydrostatic pressure ratio - ν = (c_T_inlet / K_s) ** 2 / P_0 # Tritium partial pressure ratio + ε_l = 1 - ε_g # Liquid phase fraction + ψ = (ρ_l * g * ε_l * L) / P_0 # Hydrostatic pressure ratio (Eq. 8.3) + ν = (c_T_inlet / K_s) ** 2 / P_0 # Tritium partial pressure ratio (Eq. 8.5) + Bo_l = u_l * L / (ε_l * E_l) # Bodenstein number, liquid phase (Eq. 8.9) + phi_l = a * h_l * L / u_l # Transfer units parameter, liquid phase (Eq. 8.11) + Bo_g = u_g0 * L / (ε_g * E_g) # Bodenstein number, gas phase (Eq. 8.10) + phi_g = ( + 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) + ) # Transfer units parameter, gas phase (Eq. 8.12) + + y_T2_in = P_T2_in / P_0 # Inlet tritium molar fraction in gas phase dimensionless_params = { "Bo_l": Bo_l, - "phi_l": Φ_l, + "phi_l": phi_l, "Bo_g": Bo_g, - "phi_g": Φ_g, + "phi_g": phi_g, "psi": ψ, "nu": ν, } @@ -568,23 +572,56 @@ def boundary_conditions(Sa, Sb): c_T_outlet = x_T_outlet_dimless * c_T_inlet # Gas partial pressure at outlet (xi=1) - P_outlet_gas = P_0 * (1 - dimensionless_params["psi"]) - P_T2_outlet_gas = y_T2_outlet_gas * P_outlet_gas + P_outlet = P_0 * ( + 1 - dimensionless_params["psi"] + ) # Derived from Eq. 8.4 at xi=1 + P_T2_out = y_T2_outlet_gas * P_outlet + + # Mass transfer consistency check + N_A = const.N_A # Avogadro's number, 1/mol + # Tritium molar flow rate into the column via liquid + n_T_in_liquid = c_T_inlet * Q_l * N_A # Triton/s + + # Tritium molar flow rate out of the column via liquid + n_T_out_liquid = c_T_outlet * Q_l * N_A # Tritons/s + + # Tritium molar flow rate into the column via gas + n_T2_in_gas = (P_T2_in * Q_g / (R * T)) * N_A # T2/s + n_T_in_gas = n_T2_in_gas * 2 # Triton/s + + # Calculate outlet gas volumetric flow rate (gas expands as pressure drops) + Q_g_out = (P_0 * Q_g) / P_outlet + # Tritium molar flow rate out of the column via gas + n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A # T2/s + n_T_out_gas = n_T2_out_gas * 2 # Triton/s + + T_in = n_T_in_liquid + n_T_in_gas + T_out = n_T_out_liquid + n_T_out_gas + + print( + f"Tritum in (liquid phase): {n_T_in_liquid:.4e} Tritons/s" + f", Tritium in (gas phase): {n_T_in_gas:.4e} Tritons/s" + ) - # Tritium extraction rate (mol/s) - extraction_rate = Q_l * (c_T_inlet - c_T_outlet) # mol/s + print( + f"Tritium out (liquid phase): {n_T_out_liquid:.4e} Tritons/s" + f", Tritium out (gas phase): {n_T_out_gas:.4e} Tritons/s" + ) - # Gas velocity at outlet (xi=1) - u_g_outlet = u_g0 / (1 - dimensionless_params["psi"]) + print( + f"Total Tritium in: {T_in:.4e} Tritons/s" + f", Total Tritium out: {T_out:.4e} Tritons/s" + ) results = { "extraction_efficiency [%]": efficiency * 100, - "extraction_rate [mol/s]": extraction_rate, + "c_T_inlet [mol/m^3]": c_T_inlet, "c_T_outlet [mol/m^3]": c_T_outlet, "liquid_vol_flow [m^3/s]": Q_l, - "total_gas_P_outlet [Pa]": P_outlet_gas, - "P_T2_outlet_gas [Pa]": P_T2_outlet_gas, - "gas_vol_flow_outlet [m^3/s]": u_g_outlet * A_g, + "P_T2_inlet_gas [Pa]": P_T2_in, + "P_T2_outlet_gas [Pa]": P_T2_out, + "total_gas_P_outlet [Pa]": P_outlet, + "gas_vol_flow [m^3/s]": Q_g, } else: @@ -623,10 +660,10 @@ def __init__( } super().__init__(func=self.func) - def func(self, c_T_inlet, y_T2_inlet): + def func(self, c_T_inlet, P_T2_inlet): new_params = self.params.copy() new_params["c_T_inlet"] = c_T_inlet - new_params["y_T2_in"] = y_T2_inlet + new_params["P_T2_in"] = P_T2_inlet res = solve(new_params) From 8ee0eb3c2477d5b5ae6e5895dd1fae9bcac11981 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 28 Oct 2025 14:04:34 -0400 Subject: [PATCH 04/26] output mass flow rates instead --- src/python/custom_pathsim_blocks.py | 51 ++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index cb7ff12..9fc843b 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -410,11 +410,10 @@ def update(self, t): from scipy.integrate import solve_bvp from scipy import constants as const -R = 8.314 # J/(mol·K), Universal gas constant +R = 8.314 # J/(mol·K), Universal gas constant TODO read from const def solve(params): - def solve_tritium_extraction(dimensionless_params, y_T2_in): """ Solves the BVP for tritium extraction in a bubble column. @@ -528,8 +527,6 @@ def boundary_conditions(Sa, Sb): g = params["g"] T = params["T"] - R = params.get("R", R) # Ideal gas constant, J/(mol.K) - # Calculate the superficial flow velocities A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column u_g0 = Q_g / A # m/s, superficial gas inlet velocity @@ -597,20 +594,23 @@ def boundary_conditions(Sa, Sb): T_in = n_T_in_liquid + n_T_in_gas T_out = n_T_out_liquid + n_T_out_gas + print("\n") + print(f"Solver c_T_outlet: {c_T_outlet:.4e} mol/m^3") + print(f"Solver P_T2_outlet: {P_T2_out:.4e} Pa") print( - f"Tritum in (liquid phase): {n_T_in_liquid:.4e} Tritons/s" - f", Tritium in (gas phase): {n_T_in_gas:.4e} Tritons/s" + f"Tritum in (liquid phase): {n_T_in_liquid / N_A:.4e} mol Tritons/s" + f", Tritium in (gas phase): {n_T_in_gas / N_A:.4e} mol Tritons/s" ) print( - f"Tritium out (liquid phase): {n_T_out_liquid:.4e} Tritons/s" - f", Tritium out (gas phase): {n_T_out_gas:.4e} Tritons/s" + f"Tritium out (liquid phase): {n_T_out_liquid / N_A:.4e} mol Tritons/s" + f", Tritium out (gas phase): {n_T_out_gas / N_A:.4e} mol Tritons/s" ) print( - f"Total Tritium in: {T_in:.4e} Tritons/s" - f", Total Tritium out: {T_out:.4e} Tritons/s" + f"Total Tritium in: {T_in / N_A:.4e} mol Tritons/s" + f", Total Tritium out: {T_out / N_A:.4e} mol Tritons/s" ) results = { @@ -622,6 +622,8 @@ def boundary_conditions(Sa, Sb): "P_T2_outlet_gas [Pa]": P_T2_out, "total_gas_P_outlet [Pa]": P_outlet, "gas_vol_flow [m^3/s]": Q_g, + "tritium_out_liquid [mol/s]": n_T_out_liquid / N_A, + "tritium_out_gas [mol/s]": n_T_out_gas / N_A, } else: @@ -639,7 +641,24 @@ class GLC(Function): """ def __init__( - self, P_0, L, u_l, u_g0, ε_g, ε_l, E_g, E_l, a, h_l, ρ_l, K_s, D, T, g=9.81 + self, + P_0, + L, + u_l, + u_g0, + ε_g, + ε_l, + E_g, + E_l, + a, + h_l, + ρ_l, + K_s, + Q_l, + Q_g, + D, + T, + g=9.81, ): self.params = { "P_0": P_0, @@ -654,6 +673,8 @@ def __init__( "h_l": h_l, "ρ_l": ρ_l, "K_s": K_s, + "Q_l": Q_l, + "Q_g": Q_g, "g": g, "D": D, "T": T, @@ -670,4 +691,10 @@ def func(self, c_T_inlet, P_T2_inlet): c_T_outlet = res["c_T_outlet [mol/m^3]"] P_T2_outlet = res["P_T2_outlet_gas [Pa]"] - return c_T_outlet, P_T2_outlet + n_T_out_liquid = res["tritium_out_liquid [mol/s]"] + n_T_out_gas = res["tritium_out_gas [mol/s]"] + print( + f"GLC block: c_T_outlet = {c_T_outlet:.4e} mol/m^3, P_T2_outlet = {P_T2_outlet:.4e} Pa" + ) + + return n_T_out_liquid, n_T_out_gas From 3fa39db83bebd132ae38f5c05597e2d30bbd6a9c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Oct 2025 14:00:39 -0400 Subject: [PATCH 05/26] removed print statements + raise error --- src/python/custom_pathsim_blocks.py | 23 +---------------------- 1 file changed, 1 insertion(+), 22 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 9fc843b..79c3e26 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -592,27 +592,6 @@ def boundary_conditions(Sa, Sb): n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A # T2/s n_T_out_gas = n_T2_out_gas * 2 # Triton/s - T_in = n_T_in_liquid + n_T_in_gas - T_out = n_T_out_liquid + n_T_out_gas - print("\n") - print(f"Solver c_T_outlet: {c_T_outlet:.4e} mol/m^3") - print(f"Solver P_T2_outlet: {P_T2_out:.4e} Pa") - - print( - f"Tritum in (liquid phase): {n_T_in_liquid / N_A:.4e} mol Tritons/s" - f", Tritium in (gas phase): {n_T_in_gas / N_A:.4e} mol Tritons/s" - ) - - print( - f"Tritium out (liquid phase): {n_T_out_liquid / N_A:.4e} mol Tritons/s" - f", Tritium out (gas phase): {n_T_out_gas / N_A:.4e} mol Tritons/s" - ) - - print( - f"Total Tritium in: {T_in / N_A:.4e} mol Tritons/s" - f", Total Tritium out: {T_out / N_A:.4e} mol Tritons/s" - ) - results = { "extraction_efficiency [%]": efficiency * 100, "c_T_inlet [mol/m^3]": c_T_inlet, @@ -627,7 +606,7 @@ def boundary_conditions(Sa, Sb): } else: - print("BVP solver failed to converge.") + raise RuntimeError("BVP solver did not converge.") return results From b836a8f636e58c340876a26977130de3cb018734 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 29 Oct 2025 14:05:06 -0400 Subject: [PATCH 06/26] docstrings --- src/python/custom_pathsim_blocks.py | 39 ++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 79c3e26..f0b0dc2 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -414,12 +414,12 @@ def update(self, t): def solve(params): - def solve_tritium_extraction(dimensionless_params, y_T2_in): + def solve_tritium_extraction(dimensionless_params: dict, y_T2_in: float) -> dict: """ Solves the BVP for tritium extraction in a bubble column. Args: - params (dict): A dictionary containing the dimensionless parameters: + dimensionless_params (dict): A dictionary containing the dimensionless parameters: Bo_l, phi_l, Bo_g, phi_g, psi, nu. y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-). @@ -604,19 +604,37 @@ def boundary_conditions(Sa, Sb): "tritium_out_liquid [mol/s]": n_T_out_liquid / N_A, "tritium_out_gas [mol/s]": n_T_out_gas / N_A, } - + return results else: raise RuntimeError("BVP solver did not converge.") - return results - -from pathsim.blocks import Function +class GLC(pathsim.blocks.Function): + """ + Gas Liquid Contactor model block. Inherits from Function block. + Inputs: c_T_inlet [mol/m^3], P_T2_inlet [Pa] + Outputs: n_T_out_liquid [mol/s], n_T_out_gas [mol/s] + More details about the model can be found in: https://doi.org/10.13182/FST95-A30485 -class GLC(Function): - """ - This is the docs + Args: + P_0: Operating pressure [Pa] + L: Column height [m] + u_l: Superficial liquid velocity [m/s] + u_g0: Superficial gas inlet velocity [m/s] + ε_g: Gas phase fraction [-] + ε_l: Liquid phase fraction [-] + E_g: Gas phase axial dispersion coefficient [m^2/s] + E_l: Liquid phase axial dispersion coefficient [m^2/s] + a: Specific interfacial area [m^2/m^3] + h_l: Liquid phase volumetric mass transfer coefficient [m/s] + ρ_l: Liquid density [kg/m^3] + K_s: Tritium solubility constant [mol/(m^3·Pa)^0.5] + Q_l: Liquid volumetric flow rate [m^3/s] + Q_g: Gas volumetric flow rate [m^3/s] + D: Column diameter [m] + T: Temperature [K] + g: Gravitational acceleration [m/s^2], default is 9.81 """ def __init__( @@ -672,8 +690,5 @@ def func(self, c_T_inlet, P_T2_inlet): n_T_out_liquid = res["tritium_out_liquid [mol/s]"] n_T_out_gas = res["tritium_out_gas [mol/s]"] - print( - f"GLC block: c_T_outlet = {c_T_outlet:.4e} mol/m^3, P_T2_outlet = {P_T2_outlet:.4e} Pa" - ) return n_T_out_liquid, n_T_out_gas From b2f69272c89a7f9c1264b8af5f32234809a44d17 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Oct 2025 08:37:58 -0400 Subject: [PATCH 07/26] custom styling --- src/components/nodes/BubblerNode.jsx | 66 ++++++++++++------------ src/components/nodes/GLCNode.jsx | 76 ++++++++++++++++++++++++++++ src/nodeConfig.js | 3 +- 3 files changed, 111 insertions(+), 34 deletions(-) create mode 100644 src/components/nodes/GLCNode.jsx diff --git a/src/components/nodes/BubblerNode.jsx b/src/components/nodes/BubblerNode.jsx index 2caa387..0b5fb74 100644 --- a/src/components/nodes/BubblerNode.jsx +++ b/src/components/nodes/BubblerNode.jsx @@ -20,29 +20,29 @@ export default function BubblerNode({ data }) { {/* Labels for sample in handles */} -
Sample in
-
soluble
-
@@ -52,52 +52,52 @@ export default function BubblerNode({ data }) { -
Vials 1
-
2
-
3
-
4
- + -
+
{data.label}
+ + +
+ c_T_inlet +
+ +
+ P_T2_in +
+ + + + + +
+ T_out_liquid +
+ +
+ T_out_gas +
+ + + +
+ ); +} \ No newline at end of file diff --git a/src/nodeConfig.js b/src/nodeConfig.js index 4dc91d9..41d2024 100644 --- a/src/nodeConfig.js +++ b/src/nodeConfig.js @@ -16,6 +16,7 @@ import BubblerNode from './components/nodes/BubblerNode'; import WallNode from './components/nodes/WallNode'; import { DynamicHandleNode } from './components/nodes/DynamicHandleNode'; import SwitchNode from './components/nodes/SwitchNode'; +import { GLCNode } from './components/nodes/GLCNode'; // Node types mapping export const nodeTypes = { @@ -64,7 +65,7 @@ export const nodeTypes = { ode: DynamicHandleNode, interface: DynamicHandleNode, switch: SwitchNode, - glc: createFunctionNode(2, 2), + glc: GLCNode, }; export const nodeMathTypes = { From cfa11e023e0cf98d27b19baeb4beff2b5753499c Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Oct 2025 08:38:10 -0400 Subject: [PATCH 08/26] added port aliases --- src/python/custom_pathsim_blocks.py | 9 +++++++++ src/python/pathsim_utils.py | 4 ++-- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index f0b0dc2..678dff7 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -637,6 +637,15 @@ class GLC(pathsim.blocks.Function): g: Gravitational acceleration [m/s^2], default is 9.81 """ + _port_map_in = { + "c_T_inlet": 0, + "P_T2_in": 1, + } + _port_map_out = { + "T_out_liquid": 0, + "T_out_gas": 1, + } + def __init__( self, P_0, diff --git a/src/python/pathsim_utils.py b/src/python/pathsim_utils.py index 9e9525c..9dd0380 100644 --- a/src/python/pathsim_utils.py +++ b/src/python/pathsim_utils.py @@ -519,7 +519,7 @@ def get_input_index(block: Block, edge: dict, block_to_input_index: dict) -> int return edge["targetHandle"] # TODO maybe we could directly use the targetHandle as a port alias for these: - if type(block) in (Function, ODE, pathsim.blocks.Switch, GLC): + if type(block) in (Function, ODE, pathsim.blocks.Switch): return int(edge["targetHandle"].replace("target-", "")) if isinstance(block, Adder): if block.operations: @@ -564,7 +564,7 @@ def get_output_index(block: Block, edge: dict) -> int: ) return output_index # TODO maybe we could directly use the targetHandle as a port alias for these: - if type(block) in (Function, ODE, GLC): + if type(block) in (Function, ODE): # Function and ODE outputs are always in order, so we can use the handle directly assert edge["sourceHandle"], edge return int(edge["sourceHandle"].replace("source-", "")) From 11281d789d6c530fe02f3416cae84b3840ba9808 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Oct 2025 09:06:00 -0400 Subject: [PATCH 09/26] R from scipy --- src/python/custom_pathsim_blocks.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 678dff7..5f028be 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -410,7 +410,7 @@ def update(self, t): from scipy.integrate import solve_bvp from scipy import constants as const -R = 8.314 # J/(mol·K), Universal gas constant TODO read from const +# R = 8.314 # J/(mol·K), Universal gas constant TODO read from const def solve(params): @@ -540,7 +540,7 @@ def boundary_conditions(Sa, Sb): phi_l = a * h_l * L / u_l # Transfer units parameter, liquid phase (Eq. 8.11) Bo_g = u_g0 * L / (ε_g * E_g) # Bodenstein number, gas phase (Eq. 8.10) phi_g = ( - 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) + 0.5 * (const.R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) ) # Transfer units parameter, gas phase (Eq. 8.12) y_T2_in = P_T2_in / P_0 # Inlet tritium molar fraction in gas phase @@ -583,13 +583,13 @@ def boundary_conditions(Sa, Sb): n_T_out_liquid = c_T_outlet * Q_l * N_A # Tritons/s # Tritium molar flow rate into the column via gas - n_T2_in_gas = (P_T2_in * Q_g / (R * T)) * N_A # T2/s + n_T2_in_gas = (P_T2_in * Q_g / (const.R * T)) * N_A # T2/s n_T_in_gas = n_T2_in_gas * 2 # Triton/s # Calculate outlet gas volumetric flow rate (gas expands as pressure drops) Q_g_out = (P_0 * Q_g) / P_outlet # Tritium molar flow rate out of the column via gas - n_T2_out_gas = (P_T2_out * Q_g_out / (R * T)) * N_A # T2/s + n_T2_out_gas = (P_T2_out * Q_g_out / (const.R * T)) * N_A # T2/s n_T_out_gas = n_T2_out_gas * 2 # Triton/s results = { From 38752501f3d7435e53f8c2a1c8a2ec8d6da9acac Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Oct 2025 14:27:13 -0400 Subject: [PATCH 10/26] efficiency as dynamic output + fixed negative pressure --- src/components/nodes/GLCNode.jsx | 16 ++++-- src/python/custom_pathsim_blocks.py | 75 +++++++++++++++++++---------- 2 files changed, 63 insertions(+), 28 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index 6e0c670..d1e5419 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -40,13 +40,22 @@ export function GLCNode({ data }) { fontWeight: 'bold', textAlign: 'right' }}> - P_T2_in + y_T2_in
- - + +
+ Efficiency +
+
); } \ No newline at end of file diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 5f028be..2b5044c 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -410,11 +410,9 @@ def update(self, t): from scipy.integrate import solve_bvp from scipy import constants as const -# R = 8.314 # J/(mol·K), Universal gas constant TODO read from const - def solve(params): - def solve_tritium_extraction(dimensionless_params: dict, y_T2_in: float) -> dict: + def solve_tritium_extraction(dimensionless_params, y_T2_in, elements): """ Solves the BVP for tritium extraction in a bubble column. @@ -422,6 +420,7 @@ def solve_tritium_extraction(dimensionless_params: dict, y_T2_in: float) -> dict dimensionless_params (dict): A dictionary containing the dimensionless parameters: Bo_l, phi_l, Bo_g, phi_g, psi, nu. y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-). + elements (int): Number of mesh elements for the solver. Returns: sol: The solution object from scipy.integrate.solve_bvp. @@ -451,7 +450,7 @@ def ode_system(xi, S): dS0_dxi = dx_T_dxi # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2 - dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) + dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) # Eq. 9.1.4 # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi) dS2_dxi = dy_T2_dxi @@ -461,9 +460,9 @@ def ode_system(xi, S): denominator = 1 - psi * xi denominator = np.where(np.isclose(denominator, 0), 1e-9, denominator) - term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi # Part of Eq. 9.3.3 (fourth line) - term2 = phi_g * theta # Part of Eq. 9.3.3 (fourth line) - dS3_dxi = (Bo_g / denominator) * (term1 - term2) # Eq. 9.3.3 (fourth line) + term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi + term2 = phi_g * theta + dS3_dxi = (Bo_g / denominator) * (term1 - term2) return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi)) @@ -483,7 +482,7 @@ def boundary_conditions(Sa, Sb): res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1]) # Eq. 10.2 # At xi = 0: y_T2(0) = y_T2(0-) + (1/Bo_g) * dy_T2/d(xi)|_0 - res3 = Sa[2] - (y_T2_in + (1 / Bo_g) * Sa[3]) # Eq. 10.3 + res3 = Sa[2] - y_T2_in - (1 / Bo_g) * Sa[3] # Eq. 10.3 # At xi = 1: dy_T2/d(xi) = 0 res4 = Sb[3] # Eq. 10.4 @@ -491,22 +490,32 @@ def boundary_conditions(Sa, Sb): return np.array([res1, res2, res3, res4]) # Set up the mesh and an initial guess for the solver. - xi = np.linspace(0, 1, 51) + xi = np.linspace(0, 1, elements + 1) - # A flat initial guess is often good enough to get the solver started. + # An initial guess that is physically plausible can significantly help convergence. + # We expect liquid concentration (x_T) to decrease from inlet (xi=1) to outlet (xi=0). + # We expect gas concentration (y_T2) to increase from inlet (xi=0) to outlet (xi=1). y_guess = np.zeros((4, xi.size)) - y_guess[0, :] = 0.8 # Guess for x_T - y_guess[2, :] = 1e-5 # Guess for y_T2 + y_guess[0, :] = np.linspace( + 0.5, 1.0, xi.size + ) # Guess for x_T (linear decrease) + y_guess[1, :] = -0.5 # Guess for dx_T/dxi + y_guess[2, :] = np.linspace( + y_T2_in, y_T2_in + 1e-4, xi.size + ) # Guess for y_T2 (linear increase) + y_guess[3, :] = 1e-4 # Guess for dy_T2/dxi # Run the BVP solver - sol = solve_bvp(ode_system, boundary_conditions, xi, y_guess, tol=1e-5) + sol = solve_bvp( + ode_system, boundary_conditions, xi, y_guess, tol=1e-5, max_nodes=10000 + ) return sol # Unpack parameters c_T_inlet = params["c_T_inlet"] - P_T2_in = params["P_T2_in"] - P_0 = params["P_0"] + y_T2_in = params["y_T2_in"] + P_outlet = params["P_outlet"] ρ_l = params["ρ_l"] K_s = params["K_s"] @@ -527,6 +536,16 @@ def boundary_conditions(Sa, Sb): g = params["g"] T = params["T"] + elements = params["elements"] # Number of mesh elements for solver + + # Calculate inlet pressure hydrostatically + P_0 = P_outlet + ρ_l * g * L + + if P_0 <= 0: + raise ValueError( + f"Calculated inlet pressure P_0 must be positive, but got {P_0:.2e} Pa. Check P_outlet, rho_l, g, and L." + ) + # Calculate the superficial flow velocities A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column u_g0 = Q_g / A # m/s, superficial gas inlet velocity @@ -543,8 +562,6 @@ def boundary_conditions(Sa, Sb): 0.5 * (const.R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) ) # Transfer units parameter, gas phase (Eq. 8.12) - y_T2_in = P_T2_in / P_0 # Inlet tritium molar fraction in gas phase - dimensionless_params = { "Bo_l": Bo_l, "phi_l": phi_l, @@ -555,7 +572,7 @@ def boundary_conditions(Sa, Sb): } # Solve the model - solution = solve_tritium_extraction(dimensionless_params, y_T2_in) + solution = solve_tritium_extraction(dimensionless_params, y_T2_in, elements) # --- Results --- if solution.success: @@ -583,6 +600,7 @@ def boundary_conditions(Sa, Sb): n_T_out_liquid = c_T_outlet * Q_l * N_A # Tritons/s # Tritium molar flow rate into the column via gas + P_T2_in = y_T2_in * P_0 # [Pa] n_T2_in_gas = (P_T2_in * Q_g / (const.R * T)) * N_A # T2/s n_T_in_gas = n_T2_in_gas * 2 # Triton/s @@ -592,6 +610,9 @@ def boundary_conditions(Sa, Sb): n_T2_out_gas = (P_T2_out * Q_g_out / (const.R * T)) * N_A # T2/s n_T_out_gas = n_T2_out_gas * 2 # Triton/s + T_in = n_T_in_liquid + n_T_in_gas + T_out = n_T_out_liquid + n_T_out_gas + results = { "extraction_efficiency [%]": efficiency * 100, "c_T_inlet [mol/m^3]": c_T_inlet, @@ -618,7 +639,7 @@ class GLC(pathsim.blocks.Function): More details about the model can be found in: https://doi.org/10.13182/FST95-A30485 Args: - P_0: Operating pressure [Pa] + P_outlet: Outlet operating pressure [Pa] L: Column height [m] u_l: Superficial liquid velocity [m/s] u_g0: Superficial gas inlet velocity [m/s] @@ -639,16 +660,17 @@ class GLC(pathsim.blocks.Function): _port_map_in = { "c_T_inlet": 0, - "P_T2_in": 1, + "y_T2_in": 1, } _port_map_out = { "T_out_liquid": 0, "T_out_gas": 1, + "efficiency": 2, } def __init__( self, - P_0, + P_outlet, L, u_l, u_g0, @@ -665,9 +687,10 @@ def __init__( D, T, g=9.81, + initial_nb_of_elements=20, ): self.params = { - "P_0": P_0, + "P_outlet": P_outlet, "L": L, "u_l": u_l, "u_g0": u_g0, @@ -684,13 +707,14 @@ def __init__( "g": g, "D": D, "T": T, + "elements": initial_nb_of_elements, } super().__init__(func=self.func) - def func(self, c_T_inlet, P_T2_inlet): + def func(self, c_T_inlet, y_T2_inlet): new_params = self.params.copy() new_params["c_T_inlet"] = c_T_inlet - new_params["P_T2_in"] = P_T2_inlet + new_params["y_T2_in"] = y_T2_inlet res = solve(new_params) @@ -699,5 +723,6 @@ def func(self, c_T_inlet, P_T2_inlet): n_T_out_liquid = res["tritium_out_liquid [mol/s]"] n_T_out_gas = res["tritium_out_gas [mol/s]"] + eff = res["extraction_efficiency [%]"] - return n_T_out_liquid, n_T_out_gas + return n_T_out_liquid, n_T_out_gas, eff From 0b7ba3e207582403b377d10be2c5be4f7a608c9e Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 30 Oct 2025 14:39:15 -0400 Subject: [PATCH 11/26] zero initial guess --- src/python/custom_pathsim_blocks.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 2b5044c..867f7c4 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -495,15 +495,18 @@ def boundary_conditions(Sa, Sb): # An initial guess that is physically plausible can significantly help convergence. # We expect liquid concentration (x_T) to decrease from inlet (xi=1) to outlet (xi=0). # We expect gas concentration (y_T2) to increase from inlet (xi=0) to outlet (xi=1). + # y_guess = np.zeros((4, xi.size)) + # y_guess[0, :] = np.linspace( + # 0.5, 1.0, xi.size + # ) # Guess for x_T (linear decrease) + # y_guess[1, :] = -0.5 # Guess for dx_T/dxi + # y_guess[2, :] = np.linspace( + # y_T2_in, y_T2_in + 1e-4, xi.size + # ) # Guess for y_T2 (linear increase) + # y_guess[3, :] = 1e-4 # Guess for dy_T2/dxi + + # a zero initial guess works better for low concentrations y_guess = np.zeros((4, xi.size)) - y_guess[0, :] = np.linspace( - 0.5, 1.0, xi.size - ) # Guess for x_T (linear decrease) - y_guess[1, :] = -0.5 # Guess for dx_T/dxi - y_guess[2, :] = np.linspace( - y_T2_in, y_T2_in + 1e-4, xi.size - ) # Guess for y_T2 (linear increase) - y_guess[3, :] = 1e-4 # Guess for dy_T2/dxi # Run the BVP solver sol = solve_bvp( @@ -539,13 +542,9 @@ def boundary_conditions(Sa, Sb): elements = params["elements"] # Number of mesh elements for solver # Calculate inlet pressure hydrostatically + assert P_outlet > 0, "P_outlet must be greater than zero." P_0 = P_outlet + ρ_l * g * L - if P_0 <= 0: - raise ValueError( - f"Calculated inlet pressure P_0 must be positive, but got {P_0:.2e} Pa. Check P_outlet, rho_l, g, and L." - ) - # Calculate the superficial flow velocities A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column u_g0 = Q_g / A # m/s, superficial gas inlet velocity From c7119f9a29c697ec14bf507a44c847bac8e98c38 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Tue, 4 Nov 2025 16:42:27 -0500 Subject: [PATCH 12/26] most recent model --- src/components/nodes/GLCNode.jsx | 8 +- src/python/custom_pathsim_blocks.py | 227 ++++++++++++++++------------ 2 files changed, 138 insertions(+), 97 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index d1e5419..58836bb 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -64,7 +64,7 @@ export function GLCNode({ data }) { fontWeight: 'bold', textAlign: 'right' }}> - T_out_liquid + c_T_oulet
- T_out_gas + P_T2_out_gas
- - + + ); diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 867f7c4..0ad7241 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -409,6 +409,7 @@ def update(self, t): import numpy as np from scipy.integrate import solve_bvp from scipy import constants as const +from scipy.optimize import root_scalar def solve(params): @@ -417,10 +418,9 @@ def solve_tritium_extraction(dimensionless_params, y_T2_in, elements): Solves the BVP for tritium extraction in a bubble column. Args: - dimensionless_params (dict): A dictionary containing the dimensionless parameters: + params (dict): A dictionary containing the dimensionless parameters: Bo_l, phi_l, Bo_g, phi_g, psi, nu. y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-). - elements (int): Number of mesh elements for the solver. Returns: sol: The solution object from scipy.integrate.solve_bvp. @@ -440,29 +440,32 @@ def ode_system(xi, S): S[1] = dx_T/d(xi) S[2] = y_T2 (dimensionless gas concentration) S[3] = dy_T2/d(xi) + + x_T = c_T / c_T(L+) + xi = z / L (dimensionless position) """ x_T, dx_T_dxi, y_T2, dy_T2_dxi = S - # Dimensionless driving force theta. Eq. 8.8 - theta = x_T - np.sqrt(np.maximum(0, (1 - psi * xi) * y_T2 / nu)) + # Dimensionless driving force theta. Eq. 12 + theta = x_T - np.sqrt(((1 - (psi * xi)) / nu) * y_T2) # MATCHES PAPER # Equation for d(S[0])/d(xi) = d(x_T)/d(xi) dS0_dxi = dx_T_dxi # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2 - dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) # Eq. 9.1.4 + dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi) dS2_dxi = dy_T2_dxi # Equation for d(S[3])/d(xi) = d^2(y_T2)/d(xi)^2 from eq (11) # Avoid division by zero if (1 - psi * xi) is close to zero at xi=1 - denominator = 1 - psi * xi - denominator = np.where(np.isclose(denominator, 0), 1e-9, denominator) - term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi - term2 = phi_g * theta - dS3_dxi = (Bo_g / denominator) * (term1 - term2) + term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi # Part of Eq. 9.3.3 (fourth line) + term2 = phi_g * theta # Part of Eq. 9.3.3 (fourth line) + dS3_dxi = (Bo_g / (1 - psi * xi)) * ( + term1 - term2 + ) # Eq. 9.3.3 (fourth line) return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi)) @@ -492,25 +495,11 @@ def boundary_conditions(Sa, Sb): # Set up the mesh and an initial guess for the solver. xi = np.linspace(0, 1, elements + 1) - # An initial guess that is physically plausible can significantly help convergence. - # We expect liquid concentration (x_T) to decrease from inlet (xi=1) to outlet (xi=0). - # We expect gas concentration (y_T2) to increase from inlet (xi=0) to outlet (xi=1). - # y_guess = np.zeros((4, xi.size)) - # y_guess[0, :] = np.linspace( - # 0.5, 1.0, xi.size - # ) # Guess for x_T (linear decrease) - # y_guess[1, :] = -0.5 # Guess for dx_T/dxi - # y_guess[2, :] = np.linspace( - # y_T2_in, y_T2_in + 1e-4, xi.size - # ) # Guess for y_T2 (linear increase) - # y_guess[3, :] = 1e-4 # Guess for dy_T2/dxi - - # a zero initial guess works better for low concentrations y_guess = np.zeros((4, xi.size)) # Run the BVP solver sol = solve_bvp( - ode_system, boundary_conditions, xi, y_guess, tol=1e-5, max_nodes=10000 + ode_system, boundary_conditions, xi, y_guess, tol=1e-6, max_nodes=10000 ) return sol @@ -518,56 +507,138 @@ def boundary_conditions(Sa, Sb): # Unpack parameters c_T_inlet = params["c_T_inlet"] y_T2_in = params["y_T2_in"] - P_outlet = params["P_outlet"] - ρ_l = params["ρ_l"] - K_s = params["K_s"] + P_0 = params["P_0"] L = params["L"] D = params["D"] - ε_g = params["ε_g"] - - Q_l = params["Q_l"] - Q_g = params["Q_g"] - - a = params["a"] - E_g = params["E_g"] - E_l = params["E_l"] - - h_l = params["h_l"] + flow_l = params["flow_l"] + u_g0 = params["u_g0"] g = params["g"] T = params["T"] elements = params["elements"] # Number of mesh elements for solver - # Calculate inlet pressure hydrostatically - assert P_outlet > 0, "P_outlet must be greater than zero." - P_0 = P_outlet + ρ_l * g * L + # --- Constants --- + g = const.g # m/s^2, Gravitational acceleration + R = const.R # J/(mol·K), Universal gas constant + N_A = const.N_A # 1/mol, Avogadro's number + M_LiPb = 2.875e-25 # Kg/molecule, Lipb molecular mass + + # Calculate empirical correlations + ρ_l = 10.45e3 * (1 - 1.61e-4 * T) # kg/m^3, Liquid (LiPb) density + σ_l = 0.52 - 0.11e-3 * T # N/m, Surface tension, liquid (LiPb) - gas (He) interface + μ_l = 1.87e-4 * np.exp(11640 / (R * T)) # Pa.s, Dynamic viscosity of liquid LiPb + ν_l = μ_l / ρ_l # m^2/s, Kinematic viscosity of liquid LiPb + D_T = 2.5e-7 * np.exp( + -27000 / (R * T) + ) # m^2/s, Tritium diffusion coefficient in liquid LiPb + + K_s = 2.32e-8 * np.exp( + -1350 / (R * T) + ) # atfrac*Pa^0.5, Sievert's constant for tritium in liquid LiPb + K_s = K_s * (ρ_l / (M_LiPb * N_A)) # mol/(m^3·Pa^0.5) - # Calculate the superficial flow velocities A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column - u_g0 = Q_g / A # m/s, superficial gas inlet velocity + + # Calculate the volumetric flow rates + Q_l = flow_l / ρ_l # m^3/s, Volumetric flow rate of liquid phase + Q_g = u_g0 * A # m^3/s, Volumetric flow rate of gas phase + + # Calculate the superficial flow velocities + # u_g0 = Q_g / A # m/s, superficial gas inlet velocity u_l = Q_l / A # m/s, superficial liquid inlet velocity - # Calculate dimensionless values + # Calculate Bond, Galilei, Schmidt and Froude numbers + Bn = (g * D**2 * ρ_l) / σ_l # Bond number + Ga = (g * D**3) / ν_l**2 # Galilei number + Sc = ν_l / D_T # Schmidt number + Fr = u_g0 / (g * D) ** 0.5 # Froude number + + # Calculate dispersion coefficients + E_l = (D * u_g0) / ( + (13 * Fr) / (1 + 6.5 * (Fr**0.8)) + ) # m^2/s, Effective axial dispersion coefficient, liquid phase + E_g = ( + 0.2 * D**2 + ) * u_g0 # m^2/s, Effective axial dispersion coefficient, gas phase + + # Calculate gas hold-up (phase fraction) & mass transfer coefficient + C = 0.2 * (Bn ** (1 / 8)) * (Ga ** (1 / 12)) * Fr # C = ε_g / (1 - ε_g)^4 + + def solveEqn(ε_g, C): + # Define the equation to solve + eqn = ε_g / (1 - ε_g) ** 4 - C + return eqn + + ε_g_initial_guess = 0.1 + try: + # bracket=[0.0001, 0.9999] tells it to *only* look in this range + sol = root_scalar(solveEqn, args=(C,), bracket=[0.00001, 0.99999]) + + # print(f"--- Using root_scalar (robust method) ---") + # print(f"C value was: {C}") + # if sol.converged: + # print(f"Solved gas hold-up (εg): {sol.root:.6f}") + # # Verify it + # verification = sol.root / (1 - sol.root)**4 + # print(f"Verification (should equal C): {verification:.6f}") + # else: + # print("Solver did not converge.") + + except ValueError as e: + print( + f"Solver failed. This can happen if C is so large that no solution exists between 0 and 1." + ) + print(f"Error: {e}") + + ε_g = sol.root # Gas phase fraction ε_l = 1 - ε_g # Liquid phase fraction - ψ = (ρ_l * g * ε_l * L) / P_0 # Hydrostatic pressure ratio (Eq. 8.3) - ν = (c_T_inlet / K_s) ** 2 / P_0 # Tritium partial pressure ratio (Eq. 8.5) - Bo_l = u_l * L / (ε_l * E_l) # Bodenstein number, liquid phase (Eq. 8.9) - phi_l = a * h_l * L / u_l # Transfer units parameter, liquid phase (Eq. 8.11) - Bo_g = u_g0 * L / (ε_g * E_g) # Bodenstein number, gas phase (Eq. 8.10) - phi_g = ( - 0.5 * (const.R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) - ) # Transfer units parameter, gas phase (Eq. 8.12) + + # Calculate outlet pressure hydrostatically & check non-negative + P_outlet = P_0 - (ρ_l * (1 - ε_g) * g * L) + + if P_outlet <= 0: + raise ValueError( + f"Calculated gas outlet pressure P_outlet must be positive, but got {P_outlet:.2e} Pa. Check P_0, rho_l, g, and L are realistic." + ) + + # Calculate interfacial area + d_b = ( + 26 * (Bn**-0.5) * (Ga**-0.12) * (Fr**-0.12) + ) * D # m, Mean bubble diameter AGREES WITH PAPER + a = 6 * ε_g / d_b # m^-1, Specific interfacial area, assuming spherical bubbles + + # Calculate volumetric mass transfer coefficient, liquid-gas + h_l_a = ( + D_T * (0.6 * Sc**0.5 * Bn**0.62 * Ga**0.31 * ε_g**1.1) / (D**2) + ) # Volumetric mass transfer coefficient, liquid-gas + + h_l = h_l_a / a # Mass transfer coefficient + + # Calculate dimensionless values + + # Hydrostatic pressure ratio (Eq. 8.3) # MATCHES PAPER + psi = (ρ_l * g * (1 - ε_g) * L) / P_0 + # Tritium partial pressure ratio (Eq. 8.5) # MATCHES PAPER + nu = ((c_T_inlet / K_s) ** 2) / P_0 + # Bodenstein number, liquid phase (Eq. 8.9) # MATCHES PAPER + Bo_l = u_l * L / (ε_l * E_l) + # Transfer units parameter, liquid phase (Eq. 8.11) # MATCHES PAPER + phi_l = a * h_l * L / u_l + # Bodenstein number, gas phase (Eq. 8.10) # MATCHES PAPER ASSUMING u_g0 + Bo_g = u_g0 * L / (ε_g * E_g) + # Transfer units parameter, gas phase (Eq. 8.12) # MATCHES PAPER + phi_g = 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) dimensionless_params = { "Bo_l": Bo_l, "phi_l": phi_l, "Bo_g": Bo_g, "phi_g": phi_g, - "psi": ψ, - "nu": ν, + "psi": psi, + "nu": nu, } # Solve the model @@ -638,20 +709,10 @@ class GLC(pathsim.blocks.Function): More details about the model can be found in: https://doi.org/10.13182/FST95-A30485 Args: - P_outlet: Outlet operating pressure [Pa] + P_0: Inlet operating pressure [Pa] L: Column height [m] - u_l: Superficial liquid velocity [m/s] u_g0: Superficial gas inlet velocity [m/s] - ε_g: Gas phase fraction [-] - ε_l: Liquid phase fraction [-] - E_g: Gas phase axial dispersion coefficient [m^2/s] - E_l: Liquid phase axial dispersion coefficient [m^2/s] - a: Specific interfacial area [m^2/m^3] - h_l: Liquid phase volumetric mass transfer coefficient [m/s] - ρ_l: Liquid density [kg/m^3] - K_s: Tritium solubility constant [mol/(m^3·Pa)^0.5] Q_l: Liquid volumetric flow rate [m^3/s] - Q_g: Gas volumetric flow rate [m^3/s] D: Column diameter [m] T: Temperature [K] g: Gravitational acceleration [m/s^2], default is 9.81 @@ -662,47 +723,27 @@ class GLC(pathsim.blocks.Function): "y_T2_in": 1, } _port_map_out = { - "T_out_liquid": 0, - "T_out_gas": 1, + "c_T_outlet": 0, + "P_T2_out_gas": 1, "efficiency": 2, } def __init__( self, - P_outlet, + P_0, L, - u_l, u_g0, - ε_g, - ε_l, - E_g, - E_l, - a, - h_l, - ρ_l, - K_s, - Q_l, - Q_g, + flow_l, D, T, - g=9.81, + g=const.g, initial_nb_of_elements=20, ): self.params = { - "P_outlet": P_outlet, + "P_0": P_0, "L": L, - "u_l": u_l, "u_g0": u_g0, - "ε_g": ε_g, - "ε_l": ε_l, - "E_g": E_g, - "E_l": E_l, - "a": a, - "h_l": h_l, - "ρ_l": ρ_l, - "K_s": K_s, - "Q_l": Q_l, - "Q_g": Q_g, + "flow_l": flow_l, "g": g, "D": D, "T": T, @@ -724,4 +765,4 @@ def func(self, c_T_inlet, y_T2_inlet): n_T_out_gas = res["tritium_out_gas [mol/s]"] eff = res["extraction_efficiency [%]"] - return n_T_out_liquid, n_T_out_gas, eff + return c_T_outlet, P_T2_outlet, eff From 28a6165956114342b18d38533fa9aeeecd9972c6 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Wed, 5 Nov 2025 13:44:33 -0500 Subject: [PATCH 13/26] fix for Process --- src/python/custom_pathsim_blocks.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 0ad7241..af9d922 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -37,7 +37,10 @@ def update(self, t): outputs = [None, None] outputs[self._port_map_out["inv"]] = x outputs[self._port_map_out["mass_flow_rate"]] = mass_rate - # update the outputs + + # make sure an array of size two is returned + outputs = [float(o) for o in outputs] + self.outputs.update_from_array(outputs) From aae3becb3a03a9c9bf9823c2b523ff5b6590b36b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 13:24:19 -0500 Subject: [PATCH 14/26] removed glc --- src/python/custom_pathsim_blocks.py | 364 ---------------------------- 1 file changed, 364 deletions(-) diff --git a/src/python/custom_pathsim_blocks.py b/src/python/custom_pathsim_blocks.py index 6de0388..ed8ecbf 100644 --- a/src/python/custom_pathsim_blocks.py +++ b/src/python/custom_pathsim_blocks.py @@ -399,367 +399,3 @@ def func(self, c_0, c_L): self.outputs["flux_0"] = flux_0 self.outputs["flux_L"] = flux_L return self.outputs - - -# GLC block - -import numpy as np -from scipy.integrate import solve_bvp -from scipy import constants as const -from scipy.optimize import root_scalar - - -def solve(params): - def solve_tritium_extraction(dimensionless_params, y_T2_in, elements): - """ - Solves the BVP for tritium extraction in a bubble column. - - Args: - params (dict): A dictionary containing the dimensionless parameters: - Bo_l, phi_l, Bo_g, phi_g, psi, nu. - y_T2_in (float): Inlet tritium molar fraction in the gas phase, y_T2(0-). - - Returns: - sol: The solution object from scipy.integrate.solve_bvp. - """ - - Bo_l = dimensionless_params["Bo_l"] - phi_l = dimensionless_params["phi_l"] - Bo_g = dimensionless_params["Bo_g"] - phi_g = dimensionless_params["phi_g"] - psi = dimensionless_params["psi"] - nu = dimensionless_params["nu"] - - def ode_system(xi, S): - """ - Defines the system of 4 first-order ODEs. - S[0] = x_T (dimensionless liquid concentration) - S[1] = dx_T/d(xi) - S[2] = y_T2 (dimensionless gas concentration) - S[3] = dy_T2/d(xi) - - x_T = c_T / c_T(L+) - xi = z / L (dimensionless position) - """ - x_T, dx_T_dxi, y_T2, dy_T2_dxi = S - - # Dimensionless driving force theta. Eq. 12 - theta = x_T - np.sqrt(((1 - (psi * xi)) / nu) * y_T2) # MATCHES PAPER - - # Equation for d(S[0])/d(xi) = d(x_T)/d(xi) - dS0_dxi = dx_T_dxi - - # Equation for d(S[1])/d(xi) = d^2(x_T)/d(xi)^2 - dS1_dxi = Bo_l * (phi_l * theta - dx_T_dxi) - - # Equation for d(S[2])/d(xi) = d(y_T2)/d(xi) - dS2_dxi = dy_T2_dxi - - # Equation for d(S[3])/d(xi) = d^2(y_T2)/d(xi)^2 from eq (11) - # Avoid division by zero if (1 - psi * xi) is close to zero at xi=1 - - term1 = (1 + 2 * psi / Bo_g) * dy_T2_dxi # Part of Eq. 9.3.3 (fourth line) - term2 = phi_g * theta # Part of Eq. 9.3.3 (fourth line) - dS3_dxi = (Bo_g / (1 - psi * xi)) * ( - term1 - term2 - ) # Eq. 9.3.3 (fourth line) - - return np.vstack((dS0_dxi, dS1_dxi, dS2_dxi, dS3_dxi)) - - def boundary_conditions(Sa, Sb): - """ - Defines the boundary conditions for the BVP. - Sa: solution at xi = 0 (liquid outlet) - Sb: solution at xi = 1 (liquid inlet) - """ - # Residuals that should be zero for a valid solution. - # Based on equations (16) and (17) in the paper. - - # At xi = 0: dx_T/d(xi) = 0 - res1 = Sa[1] # Eq. 10.1 - - # At xi = 1: x_T(1) = 1 - (1/Bo_l) * dx_T/d(xi)|_1 - res2 = Sb[0] - (1 - (1 / Bo_l) * Sb[1]) # Eq. 10.2 - - # At xi = 0: y_T2(0) = y_T2(0-) + (1/Bo_g) * dy_T2/d(xi)|_0 - res3 = Sa[2] - y_T2_in - (1 / Bo_g) * Sa[3] # Eq. 10.3 - - # At xi = 1: dy_T2/d(xi) = 0 - res4 = Sb[3] # Eq. 10.4 - - return np.array([res1, res2, res3, res4]) - - # Set up the mesh and an initial guess for the solver. - xi = np.linspace(0, 1, elements + 1) - - y_guess = np.zeros((4, xi.size)) - - # Run the BVP solver - sol = solve_bvp( - ode_system, boundary_conditions, xi, y_guess, tol=1e-6, max_nodes=10000 - ) - - return sol - - # Unpack parameters - c_T_inlet = params["c_T_inlet"] - y_T2_in = params["y_T2_in"] - P_0 = params["P_0"] - - L = params["L"] - D = params["D"] - - flow_l = params["flow_l"] - u_g0 = params["u_g0"] - - g = params["g"] - T = params["T"] - - elements = params["elements"] # Number of mesh elements for solver - - # --- Constants --- - g = const.g # m/s^2, Gravitational acceleration - R = const.R # J/(mol·K), Universal gas constant - N_A = const.N_A # 1/mol, Avogadro's number - M_LiPb = 2.875e-25 # Kg/molecule, Lipb molecular mass - - # Calculate empirical correlations - ρ_l = 10.45e3 * (1 - 1.61e-4 * T) # kg/m^3, Liquid (LiPb) density - σ_l = 0.52 - 0.11e-3 * T # N/m, Surface tension, liquid (LiPb) - gas (He) interface - μ_l = 1.87e-4 * np.exp(11640 / (R * T)) # Pa.s, Dynamic viscosity of liquid LiPb - ν_l = μ_l / ρ_l # m^2/s, Kinematic viscosity of liquid LiPb - D_T = 2.5e-7 * np.exp( - -27000 / (R * T) - ) # m^2/s, Tritium diffusion coefficient in liquid LiPb - - K_s = 2.32e-8 * np.exp( - -1350 / (R * T) - ) # atfrac*Pa^0.5, Sievert's constant for tritium in liquid LiPb - K_s = K_s * (ρ_l / (M_LiPb * N_A)) # mol/(m^3·Pa^0.5) - - A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column - - # Calculate the volumetric flow rates - Q_l = flow_l / ρ_l # m^3/s, Volumetric flow rate of liquid phase - Q_g = u_g0 * A # m^3/s, Volumetric flow rate of gas phase - - # Calculate the superficial flow velocities - # u_g0 = Q_g / A # m/s, superficial gas inlet velocity - u_l = Q_l / A # m/s, superficial liquid inlet velocity - - # Calculate Bond, Galilei, Schmidt and Froude numbers - Bn = (g * D**2 * ρ_l) / σ_l # Bond number - Ga = (g * D**3) / ν_l**2 # Galilei number - Sc = ν_l / D_T # Schmidt number - Fr = u_g0 / (g * D) ** 0.5 # Froude number - - # Calculate dispersion coefficients - E_l = (D * u_g0) / ( - (13 * Fr) / (1 + 6.5 * (Fr**0.8)) - ) # m^2/s, Effective axial dispersion coefficient, liquid phase - E_g = ( - 0.2 * D**2 - ) * u_g0 # m^2/s, Effective axial dispersion coefficient, gas phase - - # Calculate gas hold-up (phase fraction) & mass transfer coefficient - C = 0.2 * (Bn ** (1 / 8)) * (Ga ** (1 / 12)) * Fr # C = ε_g / (1 - ε_g)^4 - - def solveEqn(ε_g, C): - # Define the equation to solve - eqn = ε_g / (1 - ε_g) ** 4 - C - return eqn - - ε_g_initial_guess = 0.1 - try: - # bracket=[0.0001, 0.9999] tells it to *only* look in this range - sol = root_scalar(solveEqn, args=(C,), bracket=[0.00001, 0.99999]) - - # print(f"--- Using root_scalar (robust method) ---") - # print(f"C value was: {C}") - # if sol.converged: - # print(f"Solved gas hold-up (εg): {sol.root:.6f}") - # # Verify it - # verification = sol.root / (1 - sol.root)**4 - # print(f"Verification (should equal C): {verification:.6f}") - # else: - # print("Solver did not converge.") - - except ValueError as e: - print( - f"Solver failed. This can happen if C is so large that no solution exists between 0 and 1." - ) - print(f"Error: {e}") - - ε_g = sol.root # Gas phase fraction - ε_l = 1 - ε_g # Liquid phase fraction - - # Calculate outlet pressure hydrostatically & check non-negative - P_outlet = P_0 - (ρ_l * (1 - ε_g) * g * L) - - if P_outlet <= 0: - raise ValueError( - f"Calculated gas outlet pressure P_outlet must be positive, but got {P_outlet:.2e} Pa. Check P_0, rho_l, g, and L are realistic." - ) - - # Calculate interfacial area - d_b = ( - 26 * (Bn**-0.5) * (Ga**-0.12) * (Fr**-0.12) - ) * D # m, Mean bubble diameter AGREES WITH PAPER - a = 6 * ε_g / d_b # m^-1, Specific interfacial area, assuming spherical bubbles - - # Calculate volumetric mass transfer coefficient, liquid-gas - h_l_a = ( - D_T * (0.6 * Sc**0.5 * Bn**0.62 * Ga**0.31 * ε_g**1.1) / (D**2) - ) # Volumetric mass transfer coefficient, liquid-gas - - h_l = h_l_a / a # Mass transfer coefficient - - # Calculate dimensionless values - - # Hydrostatic pressure ratio (Eq. 8.3) # MATCHES PAPER - psi = (ρ_l * g * (1 - ε_g) * L) / P_0 - # Tritium partial pressure ratio (Eq. 8.5) # MATCHES PAPER - nu = ((c_T_inlet / K_s) ** 2) / P_0 - # Bodenstein number, liquid phase (Eq. 8.9) # MATCHES PAPER - Bo_l = u_l * L / (ε_l * E_l) - # Transfer units parameter, liquid phase (Eq. 8.11) # MATCHES PAPER - phi_l = a * h_l * L / u_l - # Bodenstein number, gas phase (Eq. 8.10) # MATCHES PAPER ASSUMING u_g0 - Bo_g = u_g0 * L / (ε_g * E_g) - # Transfer units parameter, gas phase (Eq. 8.12) # MATCHES PAPER - phi_g = 0.5 * (R * T * c_T_inlet / P_0) * (a * h_l * L / u_g0) - - dimensionless_params = { - "Bo_l": Bo_l, - "phi_l": phi_l, - "Bo_g": Bo_g, - "phi_g": phi_g, - "psi": psi, - "nu": nu, - } - - # Solve the model - solution = solve_tritium_extraction(dimensionless_params, y_T2_in, elements) - - # --- Results --- - if solution.success: - # --- Dimensionless Results --- - x_T_outlet_dimless = solution.y[0, 0] - efficiency = 1 - x_T_outlet_dimless - y_T2_outlet_gas = solution.y[2, -1] # y_T2 at xi=1 - - # --- Dimensional Results --- - # Liquid concentration at outlet (xi=0) - c_T_outlet = x_T_outlet_dimless * c_T_inlet - - # Gas partial pressure at outlet (xi=1) - P_outlet = P_0 * ( - 1 - dimensionless_params["psi"] - ) # Derived from Eq. 8.4 at xi=1 - P_T2_out = y_T2_outlet_gas * P_outlet - - # Mass transfer consistency check - N_A = const.N_A # Avogadro's number, 1/mol - # Tritium molar flow rate into the column via liquid - n_T_in_liquid = c_T_inlet * Q_l * N_A # Triton/s - - # Tritium molar flow rate out of the column via liquid - n_T_out_liquid = c_T_outlet * Q_l * N_A # Tritons/s - - # Tritium molar flow rate into the column via gas - P_T2_in = y_T2_in * P_0 # [Pa] - n_T2_in_gas = (P_T2_in * Q_g / (const.R * T)) * N_A # T2/s - n_T_in_gas = n_T2_in_gas * 2 # Triton/s - - # Calculate outlet gas volumetric flow rate (gas expands as pressure drops) - Q_g_out = (P_0 * Q_g) / P_outlet - # Tritium molar flow rate out of the column via gas - n_T2_out_gas = (P_T2_out * Q_g_out / (const.R * T)) * N_A # T2/s - n_T_out_gas = n_T2_out_gas * 2 # Triton/s - - T_in = n_T_in_liquid + n_T_in_gas - T_out = n_T_out_liquid + n_T_out_gas - - results = { - "extraction_efficiency [%]": efficiency * 100, - "c_T_inlet [mol/m^3]": c_T_inlet, - "c_T_outlet [mol/m^3]": c_T_outlet, - "liquid_vol_flow [m^3/s]": Q_l, - "P_T2_inlet_gas [Pa]": P_T2_in, - "P_T2_outlet_gas [Pa]": P_T2_out, - "total_gas_P_outlet [Pa]": P_outlet, - "gas_vol_flow [m^3/s]": Q_g, - "tritium_out_liquid [mol/s]": n_T_out_liquid / N_A, - "tritium_out_gas [mol/s]": n_T_out_gas / N_A, - } - return results - else: - raise RuntimeError("BVP solver did not converge.") - - -class GLC(pathsim.blocks.Function): - """ - Gas Liquid Contactor model block. Inherits from Function block. - Inputs: c_T_inlet [mol/m^3], P_T2_inlet [Pa] - Outputs: n_T_out_liquid [mol/s], n_T_out_gas [mol/s] - - More details about the model can be found in: https://doi.org/10.13182/FST95-A30485 - - Args: - P_0: Inlet operating pressure [Pa] - L: Column height [m] - u_g0: Superficial gas inlet velocity [m/s] - Q_l: Liquid volumetric flow rate [m^3/s] - D: Column diameter [m] - T: Temperature [K] - g: Gravitational acceleration [m/s^2], default is 9.81 - """ - - _port_map_in = { - "c_T_inlet": 0, - "y_T2_in": 1, - } - _port_map_out = { - "c_T_outlet": 0, - "P_T2_out_gas": 1, - "efficiency": 2, - } - - def __init__( - self, - P_0, - L, - u_g0, - flow_l, - D, - T, - g=const.g, - initial_nb_of_elements=20, - ): - self.params = { - "P_0": P_0, - "L": L, - "u_g0": u_g0, - "flow_l": flow_l, - "g": g, - "D": D, - "T": T, - "elements": initial_nb_of_elements, - } - super().__init__(func=self.func) - - def func(self, c_T_inlet, y_T2_inlet): - new_params = self.params.copy() - new_params["c_T_inlet"] = c_T_inlet - new_params["y_T2_in"] = y_T2_inlet - - res = solve(new_params) - - c_T_outlet = res["c_T_outlet [mol/m^3]"] - P_T2_outlet = res["P_T2_outlet_gas [Pa]"] - - n_T_out_liquid = res["tritium_out_liquid [mol/s]"] - n_T_out_gas = res["tritium_out_gas [mol/s]"] - eff = res["extraction_efficiency [%]"] - - return c_T_outlet, P_T2_outlet, eff From afe95689362bd5664705b92245a0b32fcbe0c48d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 13:26:51 -0500 Subject: [PATCH 15/26] GLC from pathsim-chem --- glc.py | 145 ++++++++++++++++++++++++++++++++++++ pyproject.toml | 2 +- src/python/pathsim_utils.py | 3 +- 3 files changed, 147 insertions(+), 3 deletions(-) create mode 100644 glc.py diff --git a/glc.py b/glc.py new file mode 100644 index 0000000..eed17ca --- /dev/null +++ b/glc.py @@ -0,0 +1,145 @@ +import pathsim +from pathsim import Simulation, Connection +import numpy as np +import matplotlib.pyplot as plt +import pathview + +# Create global variables +import numpy as np + +c_T_inlet = 3.09e-2 # mol/m^3 (c_T(L+)), Inlet tritium concentration in liquid just before inlet +y_T2_in = 0.0 # Inlet tritium molar fraction in gas (0 = pure purge gas) +P_0 = 5e5 # Pa, Gas total pressure at inlet +ρ_l = 9000 # kg/m^3, Liquid density +K_s = 2e-6 # Tritium Sievert's constant in liquid + +L = 3.0 # m, Height of the bubble column +D = 0.5 # m, Column diameter +ε_g = 0.04 # Gas phase fraction + +Q_l = 0.01 # m^3/s, Volumetric flow rate of liquid phase +Q_g = 0.0005 # m^3/s, Volumetric flow rate of gas phase at inlet + +a = 20 # m^-1, Specific liquid-gas interfacial area + +E_g = 0.05 # m^2/s, Effective axial dispersion coefficient, gas phase +E_l = 0.01 # m^2/s, Effective axial dispersion coefficient, liquid phase + +h_l = 1e-4 # m/s, Mass transfer coefficient, tritium liquid - gas + +g = 9.81 # m/s^2, Gravitational acceleration + +# Calculated parameters +ε_l = 1 - ε_g # Liquid phase fraction + +A = np.pi * (D / 2) ** 2 # m^2, Cross-sectional area of the column +A_l = A * ε_l # m^2, Cross-sectional area of the liquid phase +A_g = A * ε_g # m^2, Cross-sectional area of the gas phase + +u_l = Q_l / A_l # m/s, superficial liquid inlet velocity (positive for downward flow) +u_g0 = Q_g / A_g # m/s, gas inlet velocity (positive for upward flow) +R = 8.314 +T = 623 +N_A = 6.0221408e23 + +# Create blocks +blocks, events = [], [] + +p_t2_inlet_5 = pathsim.blocks.sources.Constant(value=y_T2_in) +blocks.append(p_t2_inlet_5) + +blanket_6 = pathview.custom_pathsim_blocks.Process( + residence_time=2, initial_value=1, source_term=0 +) +blocks.append(blanket_6) + +storage_7 = pathview.custom_pathsim_blocks.Process(residence_time=0, source_term=0) +blocks.append(storage_7) + +to_concentration_14 = pathsim.blocks.amplifier.Amplifier(gain=Q_l**-1) +blocks.append(to_concentration_14) + +inventories_16 = pathsim.blocks.scope.Scope( + labels=["blanket (inv)", "storage (inv)", "adder 20"] +) +blocks.append(inventories_16) + +all_t_flow_rates_18 = pathsim.blocks.scope.Scope(labels=["blanket (mass_flow_rate)"]) +blocks.append(all_t_flow_rates_18) + +adder_20_20 = pathsim.blocks.adder.Adder() +blocks.append(adder_20_20) + +glc_21_21 = pathview.custom_pathsim_blocks.GLC( + P_outlet=5e5, + L=L, + u_l=u_l, + u_g0=u_g0, + ε_g=ε_g, + ε_l=ε_l, + E_g=E_g, + E_l=E_l, + a=a, + h_l=h_l, + ρ_l=ρ_l, + K_s=K_s, + Q_l=Q_l, + Q_g=Q_g, + D=D, + T=T, + initial_nb_of_elements=20, +) +blocks.append(glc_21_21) + + +# Create events + + +# Create connections + +connections = [ + Connection(blanket_6["mass_flow_rate"], to_concentration_14[0]), + Connection(blanket_6["inv"], inventories_16[0]), + Connection(storage_7["inv"], inventories_16[1]), + Connection(blanket_6["mass_flow_rate"], all_t_flow_rates_18[0]), + Connection(blanket_6["inv"], adder_20_20[0]), + Connection(storage_7["inv"], adder_20_20[1]), + Connection(adder_20_20[0], inventories_16[2]), + Connection(to_concentration_14[0], glc_21_21["c_T_inlet"]), + Connection(p_t2_inlet_5[0], glc_21_21["y_T2_in"]), + Connection(glc_21_21["T_out_gas"], storage_7[0]), + Connection(glc_21_21["T_out_liquid"], blanket_6[0]), +] + +# Create simulation +my_simulation = Simulation( + blocks, + connections, + events=events, + Solver=pathsim.solvers.SSPRK22, + dt=1, + dt_min=1e-16, + iterations_max=200, + log=True, + tolerance_fpi=1e-10, + **{}, +) + +if __name__ == "__main__": + my_simulation.run(200) + + # Optional: Plotting results + scopes = [block for block in blocks if isinstance(block, pathsim.blocks.Scope)] + fig, axs = plt.subplots( + nrows=len(scopes), sharex=True, figsize=(10, 5 * len(scopes)) + ) + for i, scope in enumerate(scopes): + plt.sca(axs[i] if len(scopes) > 1 else axs) + time, data = scope.read() + # plot the recorded data + for p, d in enumerate(data): + lb = scope.labels[p] if p < len(scope.labels) else f"port {p}" + plt.plot(time, d, label=lb) + plt.legend() + plt.xlabel("Time") + plt.show() diff --git a/pyproject.toml b/pyproject.toml index f16925b..c6b77ed 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ classifiers = [ requires-python = ">=3.8" dependencies = [ "pathsim>=0.8.2", - "pathsim-chem", + "git+https://github.com/pathsim/pathsim-chem@glc", "matplotlib>=3.7.0", "numpy>=1.24.0", "plotly>=6.0", diff --git a/src/python/pathsim_utils.py b/src/python/pathsim_utils.py index fde2764..a634b29 100644 --- a/src/python/pathsim_utils.py +++ b/src/python/pathsim_utils.py @@ -51,9 +51,8 @@ Splitter3, FestimWall, Integrator, - GLC, ) -from pathsim_chem import Bubbler4, Splitter +from pathsim_chem import Bubbler4, Splitter, GLC import inspect NAME_TO_SOLVER = { From c82a677469ce8779aaf4ea38d7913b1c9b7b7ca2 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 13:49:08 -0500 Subject: [PATCH 16/26] fixed syntax --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index c6b77ed..462ba47 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ classifiers = [ requires-python = ">=3.8" dependencies = [ "pathsim>=0.8.2", - "git+https://github.com/pathsim/pathsim-chem@glc", + "pathsim-chem@git+https://github.com/pathsim/pathsim-chem#egg=glc" "matplotlib>=3.7.0", "numpy>=1.24.0", "plotly>=6.0", From 67decf63177a0ccf32378e55a88a45803c9bdf61 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 14:00:21 -0500 Subject: [PATCH 17/26] added comma --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 462ba47..131250a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ classifiers = [ requires-python = ">=3.8" dependencies = [ "pathsim>=0.8.2", - "pathsim-chem@git+https://github.com/pathsim/pathsim-chem#egg=glc" + "pathsim-chem@git+https://github.com/pathsim/pathsim-chem#egg=glc", "matplotlib>=3.7.0", "numpy>=1.24.0", "plotly>=6.0", From a9b757dedb7bbb54beb266819ea6c2c6a8f6fa3b Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 14:01:23 -0500 Subject: [PATCH 18/26] fixed branch --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 131250a..738d537 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ classifiers = [ requires-python = ">=3.8" dependencies = [ "pathsim>=0.8.2", - "pathsim-chem@git+https://github.com/pathsim/pathsim-chem#egg=glc", + "pathsim-chem@git+https://github.com/pathsim/pathsim-chem@glc", "matplotlib>=3.7.0", "numpy>=1.24.0", "plotly>=6.0", From 7f27a67ddb4670778033d1176ffa7dbcbe51cca3 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 14:52:37 -0500 Subject: [PATCH 19/26] added new handle --- src/components/nodes/GLCNode.jsx | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index 58836bb..952e8e4 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -75,11 +75,23 @@ export function GLCNode({ data }) { fontWeight: 'bold', textAlign: 'right' }}> - P_T2_out_gas + y_T2_out + + +
+ P_out_gas
- + + ); From 0b63561b5ea29f59d25ba6f2fc4eadd51f969173 Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Thu, 13 Nov 2025 17:40:50 -0500 Subject: [PATCH 20/26] added n_T_out_gas handle --- src/components/nodes/GLCNode.jsx | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index 952e8e4..a39d917 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -4,8 +4,8 @@ export function GLCNode({ data }) { return (
+
+ n_T_out_gas +
+ +
); From 0fff964e3b070254471525a749d283b84078c44d Mon Sep 17 00:00:00 2001 From: RemDelaporteMathurin Date: Mon, 17 Nov 2025 13:28:46 -0500 Subject: [PATCH 21/26] rm print statemetn --- src/python/pathsim_utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/python/pathsim_utils.py b/src/python/pathsim_utils.py index a634b29..2cb9274 100644 --- a/src/python/pathsim_utils.py +++ b/src/python/pathsim_utils.py @@ -280,7 +280,6 @@ def make_solver_params(solver_prms: dict, eval_namespace=None): # TODO get the default from pathsim._constants prms[k] = None else: - print(v, type(v)) prms[k] = eval(v, eval_namespace) elif k == "log": if v == "true": From 2bbe572efdeb948604a2d69496aa35b0219689c0 Mon Sep 17 00:00:00 2001 From: rossmacdonald98 Date: Mon, 17 Nov 2025 21:18:13 +0000 Subject: [PATCH 22/26] Updated handles on GLC block --- src/components/nodes/GLCNode.jsx | 38 +++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index a39d917..f716c68 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -54,34 +54,34 @@ export function GLCNode({ data }) { fontWeight: 'bold', textAlign: 'right' }}> - Efficiency + c_T_out
- c_T_oulet + y_T2_out
- y_T2_out + efficiency
+ Q_l +
+ +
- n_T_out_gas + Q_g_out
- - - - - + + + + + + ); } \ No newline at end of file From aeb8f3b1af5ed669e9ce0ace9ef91a504f1f92e5 Mon Sep 17 00:00:00 2001 From: rossmacdonald98 Date: Mon, 17 Nov 2025 21:43:39 +0000 Subject: [PATCH 23/26] Updated handles on GLC block 2 --- src/components/nodes/GLCNode.jsx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index f716c68..27ebffe 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -29,7 +29,7 @@ export function GLCNode({ data }) { fontWeight: 'bold', textAlign: 'right' }}> - c_T_inlet + c_T_in
- +
- P_out_gas + P_out
- + - - + +
From 69df131aa01176d869cb6fd041b1d304ef028884 Mon Sep 17 00:00:00 2001 From: rossmacdonald98 Date: Mon, 17 Nov 2025 22:10:45 +0000 Subject: [PATCH 24/26] Updated GLC block handle positions --- src/components/nodes/GLCNode.jsx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index 27ebffe..ff56f54 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -24,7 +24,7 @@ export function GLCNode({ data }) {
Date: Tue, 18 Nov 2025 14:37:24 +0000 Subject: [PATCH 25/26] Added n_T out handles to GLC block --- src/components/nodes/GLCNode.jsx | 54 +++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index ff56f54..ac35806 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -5,7 +5,7 @@ export function GLCNode({ data }) {
- y_T2_out + efficiency
- efficiency + Q_l
- P_out + n_T_out_liquid
- Q_l + y_T2_out
+ P_out +
+ +
+
+ n_T_out_gas +
+ - - - - - + + + + + + +
); } \ No newline at end of file From f8fbd9d382d5f10c1f2dbc43c23cc8ad86b1432c Mon Sep 17 00:00:00 2001 From: Ross MacDonald Date: Thu, 20 Nov 2025 22:13:52 +0000 Subject: [PATCH 26/26] added flow handles --- src/components/nodes/GLCNode.jsx | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/src/components/nodes/GLCNode.jsx b/src/components/nodes/GLCNode.jsx index ac35806..14d2bca 100644 --- a/src/components/nodes/GLCNode.jsx +++ b/src/components/nodes/GLCNode.jsx @@ -24,7 +24,7 @@ export function GLCNode({ data }) {
+ flow_l +
+ +
- - +
+ flow_g +
+ + + + +