diff --git a/docs/_static/unit_models/BPEDdiagram.png b/docs/_static/unit_models/BPEDdiagram.png new file mode 100644 index 0000000000..39581fecf4 Binary files /dev/null and b/docs/_static/unit_models/BPEDdiagram.png differ diff --git a/docs/technical_reference/unit_models/bipolar_electrodialysis_0D.rst b/docs/technical_reference/unit_models/bipolar_electrodialysis_0D.rst new file mode 100644 index 0000000000..c759ba8bdb --- /dev/null +++ b/docs/technical_reference/unit_models/bipolar_electrodialysis_0D.rst @@ -0,0 +1,365 @@ +Bipolar Electrodialysis (0D) +============================ + +Introduction +------------ + +Bipolar electrodialysis, an electrochemical separation technology, is primarily used to generate acids and bases +from waste salts. Recently, multiple proof of concept studies have also shown that starting from Lithium Chloride solution bipolar membranes can produce Lithium Hydroxide. +These are critical for batteries. In water treatment plants starting with waste brine, at the end of water purification, yields high concentrations sodium hydroxide. +These can be new revenue streams. To produce products from salts there is a cation exchange memrbane (CEM) and a anion exchange membrane (AEM) in parallel with the bipolar membrane. +To simplify the model we do NOT account for these salt ion fluxes. This allows direct comparisons with models presented in works of Mareev et al. (2020) and Melnikov (2022). Users may choose to add ions via the feed. +In the future we will add a unit model with CEM and AEM in parallel with a bipolar membrane. A sketch of the bipolar membrane cell stack is shown in Figure 1 with **basic** and **acidic** channels, that produce base and acid +respectively. More overview of the bipolar electrodialysis technology can be found in the *References*. + +.. figure:: ../../_static/unit_models/BPEDdiagram.png + :width: 600 + :align: center + + Figure 1. Schematic representation of a bipolar electrodialysis unit + + +One bipolar membrane along with the **acidic** and **basic** channels can thus be treated as a modelling unit that can +multiply to larger-scale systems. The presented bipolar electrodialysis model establishes mathematical descriptions of +ion and water transport across the membrane along with water splitting. Modelled transfer mechanisms include +electrical migration, diffusion of ions, osmosis, electroosmosis, and water splitting. The following are the key +assumptions made: + +* The **acidic** and **basic** side channels have identical geometry. +* For each channel, component fluxes are uniform in the bulk solutions (the 0-dimensional assumption) and are set as the average of inlet and outlet of each channel. +* Steady state: all variables are independent of time. +* Co-current flow operation. +* Ideality assumptions: activity, osmotic, and van't Hoff coefficients are set at one. +* All ion-exchange membrane properties (ion and water transport number, resistance, permeability) are constant. +* Detailed concentration gradient effect at membrane-water interfaces is neglected. +* Constant pressure and temperature through each channel. +* No boundary layer, electric double layer or diffusion layer, has been considered. + +Control Volumes +--------------- + +This model has two control volumes for the acidic and basic channels. + +* **acidic** channel +* **basic** channel + +Ports +----- + +On the two control volumes, this model provides four ports (Pyomo notation in parenthesis): + +* inlet_acidic (inlet) +* outlet_acidic (outlet) +* inlet_basic (inlet) +* outlet_basic (outlet) + +Sets +---- +This model can simulate the water splitting and the transport of multiple species. All solution components +( H\ :sub:`2`\ O, neutral solutes, and ions, including Proton and Hydroxide ion) form a Pyomo set in the model. +For a clear model demonstration, this document uses a NaCl water solution along with the products of water splitting, H\ :sup:`+` and OH\ :sup:`-`, hereafter. + +This model can mathematically take a multi-component (i.e., one salt molecule to be treated) as an input; nevertheless +a multi-component solution creates unknown or difficult-to-specify parameters, e.g., the electrical transport numbers through membranes, +the multi-ion diffusivity, etc., and physical relationships, which may result in ill-posed or ill-conditioned problems challenging the models' +numerical solutions. While we continuously work on advancing our models to absorb new principles revealed by progressing +research, we advise the users be very **cautious** with simulating multi-component system by this programmed model for aspects stated above. +This unit model works with the MCAS property model. + +.. csv-table:: **Table 1.** List of Set + :header: "Description", "Symbol", "Indices" + + + "Time", ":math:`t`", "[t] ([0])\ :sup:`1`" + "Phase", ":math:`p`", "['Liq']" + "Component", ":math:`j`", "['H\ :sub:`2` \O', 'Na\ :sup:`+`', 'Cl\ :sup:`-`', 'H\ :sup:`+`', 'OH\ :sup:`-`']" + "Ion", ":math:`j`", "['Na\ :sup:`+`', 'Cl\ :sup:`-`', 'H\ :sup:`+`', 'OH\ :sup:`-`'] \ :sup:`2`" + "Membrane", "n/a", "['bpem']" + +**Notes** + :sup:`1` The time set index is set as [0] in this steady-state model and is reserved majorly for the future extension + to a dynamic model. + + :sup:`2` "Ion" is a subset of "Component" and uses the same symbol j. + + +Degrees of Freedom +------------------ +The bipolar membrane model has multiple degrees of freedom, among which temperature, pressure, and component molar flow +rate are state variables that are fixed as initial conditions. The rest are parameters that should be provided in order +to fully solve the model. The exact degrees of freedom depend on the mode of operation. For the simplest case where no water +splitting occurs and the bipolar membrane acts like a simple electrodialysis membrane these are: + +.. csv-table:: **Table 2.** List of Degree of Freedom (DOF) + :header: "Description", "Symbol", "Variable Name", "Index", "Units", "DOF Number \ :sup:`1`" + + "Temperature, inlet_acidic", ":math:`T^{acidic}`", "temperature", "None", ":math:`K`", 1 + "Temperature, inlet_basic", ":math:`T^{basic}`", "temperature", "None", ":math:`K`", 1 + "Pressure, inlet_acidic",":math:`p^{acidic}`", "temperature", "None", ":math:`Pa`", 1 + "Pressure, inlet_basic",":math:`p^{basic}`", "temperature", "None", ":math:`Pa`", 1 + "Component molar flow rate, inlet_acidic", ":math:`N_{j,in}^{acidic}`", "flow_mol_phase_comp", "[t], ['Liq'], ['H\ :sub:`2`\O', 'Na\ :sup:`+`', '\Cl\ :sup:`-`', 'H\ :sup:`+`', 'OH\ :sup:`-`']", ":math:`mol \, s^{-1}`", 5 + "Component molar flow rate, inlet_basic", ":math:`N_{j, in}^{basic}`", "flow_mol_phase_comp", "[t], ['Liq'], ['H\ :sub:`2`\O', 'Na\ :sup:`+`', '\Cl\ :sup:`-`', 'H\ :sup:`+`', 'OH\ :sup:`-`']", ":math:`mol \, s^{-1}`", 5 + "Water transport number", ":math:`t_w`", "water_trans_number_membrane", "['bpem']", "dimensionless", 1 + "Water permeability", ":math:`L`", "water_permeability_membrane", "['bpem']", ":math:`m^{-1}s^{-1}Pa^{-1}`", 1 + "Voltage or Current \ :sup:`2`", ":math:`U` or :math:`I`", "voltage or current", "[t]", ":math:`\text{V}` or :math:`A`", 1 + "Electrode areal resistance", ":math:`r_{el}`", "electrodes_resistance", "[t]", ":math:`\Omega m^2`", 1 + "Cell number", ":math:`n`", "cell_num", "None", "dimensionless", 1 + "Current utilization coefficient", ":math:`\xi`", "current_utilization", "None", "dimensionless", 1 + "Shadow factor", ":math:`\beta`", "shadow_factor", "None", "dimensionless", 1 + "Spacer thickness", ":math:`s`", "spacer_thickness", "none", ":math:`m` ", 1 + "Membrane areal resistance", ":math:`r`", "membrane_surface_resistance", "['acidic', 'basic']", ":math:`\Omega m^2`", 2 + "Cell width", ":math:`b`", "cell_width", "None", ":math:`\text{m}`", 1 + "Cell length", ":math:`l`", "cell_length", "None", ":math:`\text{m}`", 1 + "Thickness of ion exchange membranes", ":math:`\delta`", "membrane_thickness", "['bpem']", ":math:`m`", 1 + "transport number of ions in the membrane phase", ":math:`t_j`", "ion_trans_number_membrane", "['bpem'], ['Na\ :sup:`+`', '\Cl\ :sup:`-`', 'H\ :sup:`+`', 'OH\ :sup:`-`']", "dimensionless", 4 + +**Note** + :sup:`1` DOF number takes account of the indices of the corresponding parameter. + + :sup:`2` A user should provide either current or voltage as the electrical input, in correspondence to the "Constant_Current" or "Constant_Voltage" treatment mode + + +Solution component information +------------------------------ +To fully construct solution properties, users need to provide basic component information of the feed solution to use this model. Below is a sample: + +.. code-block:: + + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + +This model, by default, uses H\ :sub:`2`\ O as the solvent of the feed solution. Please note that H\ :sup:`+` and OH\ :sup:`-` information must be supplied. Otherwise an error will be thrown. + +Information regarding the property package this unit model relies on can be found here: + +:py:mod:`watertap.property_models.ion_DSPMDE_prop_pack` + +Operation without catalyst +-------------------------- + +The simplest water splitting mode is without any catalyst. Hence default the config ``has_catalyst`` is set to false. The Mass balance equations are summarized in **Table3**. Further details on these can be found in the *References*. + +.. csv-table:: **Table 3** Mass Balance Equations + :header: "Description", "Equation", "Index set" + + "Component mass balance", ":math:`N_{j, in}^{acidic \: or\: basic}-N_{j, out}^{acidic\: or\: basic}+J_j^{acidic\: or\: basic} bl=0`", ":math:`j \in \left['H_2 O', '{Na^+} ', '{Cl^-} '\right]`" + "mass transfer flux, basic, solute", ":math:`J_j = -t_j^{bpem}\frac{\xi i_{lim}}{ z_j F}`", ":math:`j \in \left['{Na_+} ', '{Cl^-} '\right]`" + "mass transfer flux, acidic, proton", ":math:`J_j = J_{diss}`", ":math:`j \in \left['{H^+} '\right]`" + "mass transfer flux, acidic, hydroxide", ":math:`J_j = 0`", ":math:`j \in \left['{OH^-} '\right]`" + "mass transfer flux, basic, proton", ":math:`J_j = 0`", ":math:`j \in \left['{H^+} '\right]`" + "mass transfer flux, basic, hydroxide", ":math:`J_j = J_{diss}`", ":math:`j \in \left['{OH^-} '\right]`" + "mass transfer flux, acidic H\ :sub:`2`\ O", ":math:`J_j = t_w^{bpem} \left(\frac{i}{F}\right)+\left(L^{bpem} \right)\left(p_{osm}^{CEM}-p_{osm}^{AEM} \right)\left(\frac{\rho_w}{M_w}\right) - 0.5 J_{diss}`", ":math:`j \in \left['H_2 O'\right]`" + "mass transfer flux, basic, H\ :sub:`2`\ O", ":math:`J_j = -t_w^{bpem} \left(\frac{i}{F}\right)-\left(L^{bpem} \right)\left(p_{osm}^{CEM}-p_{osm}^{AEM} \right)\left(\frac{\rho_w}{M_w}\right) - 0.5 J_{diss}`", ":math:`j \in \left['H_2 O'\right]`" + +Overcoming the limiting current corresponds to a potential barrier, :math:`U_{diss}`. Important quantities are either taken as user input or computed. The appropriate configurations are ``limiting_current_density_method_bpem" for limiting current or ``limiting_potential_method_bpem`` for potential barrier. +These relationships are given in **Table 4** + + + +.. csv-table:: **Table 4** Essential equations + :header: "Description", "Equation", "Condition" + + "Flux due to hydrolysis reaction", ":math:`J_{diss} = (i - i_{lim})/F`", " " + "Limiting current density", ":math:`i_{lim} =` user input constant", "``limiting_current_density_method_bpem =LimitingCurrentDensitybpemMethod.InitialValue``" + " ", ":math:`i_{lim} = D F (C_{acidic}+C_{basic})^2 / (\sigma \delta)`", "``limiting_current_density_method_bpem =LimitingCurrentDensitybpemMethod.Empirical``" + "Potential barrier",":math:`U_{diss} =` user input constant", "``limiting_potential_method_bpem =LimitingpotentialMethod.InitialValue``" + " ", ":math:`U_{diss} = E_{crit}\lambda`", "``limiting_potential_method_bpem =LimitingpotentialMethod.Empirical``" + "Depletion length", ":math:`\lambda = E_{crit} \epsilon_0 \epsilon_r / (F \sigma)`", "``limiting_potential_method_bpem =LimitingpotentialMethod.Empirical``" + "Water splitting rate at electric field :math:`E` ", ":math:`R_{H^+/OH^-} (E) = [k_2(0)f(E)C_{H_2O}-k_r C_{H^+}C_{OH^-} ]`", "``limiting_potential_method_bpem =LimitingpotentialMethod.InitialValue``" + "Critical electric field", ":math:`R_{H^+/OH^-}(E = E_{crit})F/\lambda= 0.1 i_{lim}`", "``limiting_potential_method_bpem =LimitingpotentialMethod.Empirical``" + + +The quantities :math:`C_{H_2 O}, C_{H^+}, C_{OH^-}` are the water proton and hydroxyl concentration in +:math:`mol\, m^{-3}` and are taken to be constants. :math:`f(E)` is the second Wien effect driven enhancement of the +dissociation rate under applied electric field. It requires as input temperature and relative permittivity (:math:`\epsilon_r`). +Please note that since the unit model is assumed to operate in the water splitting regime and so :math:`U_{diss}` is always computed when ``has_catalyst`` is False. +It should be noted that :math:`J_{diss}` arises from the hydrolysis reaction and can be equated to the reaction rate. +However, as Wilhelm et al. (2001) have shown, this can be simplified to the current in excess of the limiting current. + + +.. csv-table:: **Table 5** DOF for water splitting without catalyst + :header: "Description", "Symbol", "Variable Name", "Index", "Units" + + "Diffusivity", ":math:`D`", "diffus_mass", "[bpem]", ":math:`m^2 s^{-1}`" + "Salt concentration, basic side ", ":math:`C_{basic}`", "salt_conc_basic", "[bpem]",":math:`mol m^{-3}`" + "Salt concentration, acidic side ", ":math:`C_{acidic}`", "salt_conc_acidic", "[bpem]",":math:`mol m^{-3}`" + "Membrane Fixed charge ", ":math:`\sigma`", "membrane_fixed_charge", "[bpem]",":math:`mol m^{-3}`" + "Dissociation rate constant, zero electric field ", ":math:`k_2(0)`", "kd_zero", "[bpem]",":math:`s^{-1}`" + "Recombination rate constant ", ":math:`k_r`", "kr", "[bpem]",":math:`L^1 mol^{-1} s^{-1}`" + "Relative permittivity ", ":math:`\epsilon_r`", "relative_permittivity", "[bpem]","Non-dimensional" + +.. csv-table:: **Table 6** Electrical and Performance Equations + :header: "Description", "Equation" + + "Current density", ":math:`i = \frac{I}{\beta bl}`" + "Potential drop", ":math:`U = n U_{diss} + i r_{tot}`" + "Resistance calculation", ":math:`r_{tot}=n\left(r^{acidic}+r^{basic}\right)+r_{el}`" + "Electrical power consumption", ":math:`P=UI`" + +All equations are coded as "constraints" (Pyomo). Isothermal and isobaric conditions apply. + +The model used here is derived from works by Wilhelm et al. (2002) and Ionescu, Viorel (2023).It has been validated using the bipolar membrane information available online: Fumatech, Technical Data Sheet for +Fumasep FBM, 2020. Additional inputs were obtained from from Ionescu, Viorel (2023). + + +Operation with catalyst +-------------------------- + +Choosing config ``has_catalyst`` to True enables catalyst action. With catalyst present the Mass balance term still follows the equations in **Table 3**. + + +The flux from water splitting with catalyst action is shown in **Table 7** + +.. csv-table:: **Table 7** Essential equations + :header: "Description", "Equation" + + "Water splitting flux", ":math:`J_{diss} =R_{K_A} \lambda + R_{K_B} \lambda`" + "Water splitting rate", ":math:`R_{K_A/K_B} = \frac{Q_m}{K_{A/B}}[k_2(0)f(E)C_{H_2O}-k_r C_{H^+}C_{OH^-} ]`" + "Depletion length", ":math:`\lambda = E \epsilon_0 \epsilon_r / (F \sigma)`" + "Electric current density", ":math:`i = i_{lim} + F J_{diss}`" + "Potential drop", ":math:`U=n E/\lambda + i r_{tot}`" + +Please note that since the unit model is assumed to operate in the water splitting regime and so :math:`i_{lim}` is always computed when ``has_catalyst`` is True. + + +The parameters used are given in **Table 8**. + +.. csv-table:: **Table 8.** DOF for water splitting with catalyst + :header: "Description", "Symbol", "Variable Name", "Index", "Units" + + "Catalyst concentration on the cation exchange side", ":math:`Q_m`", "membrane_fixed_catalyst_cem", "[bpem]", ":math:`mol \, m^{-3}`" + "Catalyst concentration on the anion exchange side", ":math:`Q_m`", "membrane_fixed_catalyst_aem", "[bpem]", ":math:`mol \, m^{-3}`" + "Equilibrium constant of proton disassociation", ":math:`K_A`", "k_a", "none",":math:`mol \, m^{-3}`" + "Equilibrium constant of hydroxide disassociation", ":math:`K_B`", "k_b", "none",":math:`mol \, m^{-3}`" + +The model used here is based on the analysis by Mareev et al. (2020). It and has been validated using the experimental data on bipolar membrane information available in Wilhelm et al. (2002). Additionaly inputs were obtained from Mareev et al. (2020). + +Frictional pressure drop +^^^^^^^^^^^^^^^^^^^^^^^^ +This model can optionally calculate pressured drops along the flow path in the diluate and concentrate channels through +config ``has_pressure_change`` and ``pressure_drop_method``. Under the assumption of identical diluate and concentrate +channels and starting flow rates, the flow velocities in the two channels are approximated equal and invariant over the +channel length when calculating the frictional pressure drops. This approximation is based on the evaluation that the +actual velocity variation over the channel length caused by water mass transfer across the consecutive channels leads to +negligible errors as compared to the uncertainties carried by the frictional pressure method itself. **Table 9** gives +essential equations to simulate the pressure drop. Among extensive literatures using these equations, a good reference +paper is by Wright et. al., 2018 (*References*). + +.. csv-table:: **Table 9** Essential equations supporting the pressure drop calculation + :header: "Description", "Equation", "Condition" + + "Frictional pressure drop, Darcy_Weisbach", ":math:`p_L=f\frac{\rho v^2}{2d_H}` \ :sup:`1`", "`has_pressure_change == True` and `pressure_drop_method == PressureDropMethod.Darcy_Weisbach`" + " ", ":math:`p_L=` user-input constant", "`has_pressure_change == True` and `pressure_drop_method == PressureDropMethod.Experimental`" + "Hydraulic diameter", ":math:`d_H=\frac{2db(1-\epsilon)}{d+b}`", "`hydraulic_diameter_method == HydraulicDiameterMethod.conventional`" + " ", ":math:`d_H=\frac{4\epsilon}{\frac{2}{h}+(1-\epsilon)S_{v,sp}}`", "`hydraulic_diameter_method == HydraulicDiameterMethod.spacer_specific_area_known`" + "Reynold number", ":math:`Re=\frac{\rho v d_H}{\mu}`", "`has_pressure_change == True` or `limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical`" + "Schmidt number", ":math:`Sc=\frac{\mu}{\rho D_b}`", "`has_pressure_change == True` or `limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical`" + "Sherwood number", ":math:`Sh=0.29Re^{0.5}Sc^{0.33}`", "`has_pressure_change == True` or `limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical`" + "Darcy's frictional factor", ":math:`f=4\times 50.6\epsilon^{-7.06}Re^{-1}`", "`friction_factor_method == FrictionFactorMethod.Gurreri`" + " ", ":math:`f=4\times 9.6 \epsilon^{-1} Re^{-0.5}`", "`friction_factor_method == FrictionFactorMethod.Kuroda`" + "Pressure balance", ":math:`p_{in}-p_L l =p_{out}`", "`has_pressure_change == True`" + +**Note** + + :sup:`1` We assumed a constant linear velocity (in the cell length direction), :math:`v`, in both channels and along the flow path. This :math:`v` is calculated based on the average of inlet and outlet volumetric flow rate. + +Nomenclature +------------ +.. csv-table:: **Table 10** Nomenclature + :header: "Symbol", "Description", "Unit" + :widths: 10, 20, 10 + + "**Parameters**" + ":math:`\rho_w`", "Mass density of water", ":math:`kg\ m^{-3}`" + ":math:`M_w`", "Molecular weight of water", ":math:`kg\ mol^{-1}`" + "**Variables and Parameters**" + ":math:`N`", "Molar flow rate of a component", ":math:`mol\ s^{-1}`" + ":math:`J`", "Molar flux of a component", ":math:`mol\ m^{-2}s^{-1}`" + ":math:`b`", "Cell/membrane width", ":math:`m`" + ":math:`l`", "Cell/membrane length", ":math:`m`" + ":math:`t`", "Ion transport number", "dimensionless" + ":math:`I`", "Current", ":math:`A`" + ":math:`i`", "Current density", ":math:`A m^{-2}`" + ":math:`U`", "Voltage over a stack", ":math:`V`" + ":math:`n`", "Cell number", "dimensionless" + ":math:`\xi`", "Current utilization coefficient (including ion diffusion and water electroosmosis)", "dimensionless" + ":math:`\beta`", "Shadow factor", "dimensionless" + ":math:`z`", "Ion charge", "dimensionless" + ":math:`F`", "Faraday constant", ":math:`C\ mol^{-1}`" + ":math:`\epsilon_0`", "permittivity of free space", ":math:`C\ mol^{-1}`" + ":math:`D`", "Ion Diffusivity", ":math:`F m^-1`" + ":math:`\delta`", "Membrane thickness", ":math:`m`" + ":math:`c`", "Solute concentration", ":math:`mol\ m^{-3}`" + ":math:`t_w`", "Water electroosmotic transport number", "dimensionless" + ":math:`L`", "Water permeability (osmosis)", ":math:`ms^{-1}Pa^{-1}`" + ":math:`p_{osm}`", "Osmotic pressure", ":math:`Pa`" + ":math:`r_{tot}`", "Total areal resistance", ":math:`\Omega m^2`" + ":math:`r`", "Membrane areal resistance", ":math:`\Omega m^2`" + ":math:`r_{el}`", "Electrode areal resistance", ":math:`\Omega m^2`" + ":math:`d`", "Spacer thickness", ":math:`m`" + ":math:`P`", "Power consumption", ":math:`W`" + ":math:`Q`", "Volume flow rate", ":math:`m^3s^{-1}`" + ":math:`\phi_d^{ohm}`", "Ohmic potential across a Nernst diffusion layer", ":math:`V`" + "**Subscripts and superscripts**" + ":math:`j`", "Component index", + ":math:`in`", "Inlet", + ":math:`out`", "Outlet", + ":math:`acidic`", "Cation exchange side of bipolar membrane", + ":math:`basic`", "Anion exchange side of bipolar membrane", + +Class Documentation +------------------- + +* :mod:`watertap.unit_models.Bipolar_Electrodialysis_0D` + +References +---------- +Campione, A., Gurreri, L., Ciofalo, M., Micale, G., Tamburini, A., & Cipollina, A. (2018). +Electrodialysis for water desalination: A critical assessment of recent developments on process +fundamentals, models and applications. Desalination, 434, 121-160. + +Campione, A., Cipollina, A., Bogle, I. D. L., Gurreri, L., Tamburini, A., Tedesco, M., & Micale, G. (2019). +A hierarchical model for novel schemes of electrodialysis desalination. Desalination, 465, 79-93. + +Fumatech, Technical Data Sheet for Fumasep FBM, 2020. + +Ionescu, V., 2023, March. A simple one-dimensional model for analysis of a bipolar membrane used in electrodialysis desalination. In Advanced Topics in Optoelectronics, Microelectronics, and Nanotechnologies XI (Vol. 12493, pp. 520-529). SPIE. + +Mareev, S.A., Evdochenko, E., Wessling, M., Kozaderova, O.A., Niftaliev, S.I., Pismenskaya, N.D. and Nikonenko, V.V., 2020. A comprehensive mathematical model of water splitting in bipolar membranes: Impact of the spatial distribution of fixed charges and catalyst at bipolar junction. Journal of Membrane Science, 603, p.118010. + +Melnikov, S., 2022. Ion Transport and Process of Water Dissociation in Electromembrane System with Bipolar Membrane: Modelling of Symmetrical Case. Membranes, 13(1), p.47. + +Spiegler, K. S. (1971). Polarization at ion exchange membrane-solution interfaces. Desalination, 9(4), 367-385. + +Strathmann, H. (2004). Ion-exchange membrane separation processes. Elsevier. Ch. 4. + +Strathmann, H. (2010). Electrodialysis, a mature technology with a multitude of new applications. +Desalination, 264(3), 268-288. + +Wright, N. C., Shah, S. R., & Amrose, S. E. (2018). +A robust model of brackish water electrodialysis desalination with experimental comparison at different size scales. +Desalination, 443, 27-43. + +Wilhelm, F.G., Pünt, I., Van Der Vegt, N.F.A., Wessling, M. and Strathmann, H., 2001. Optimisation strategies for the preparation of bipolar membranes with reduced salt ion leakage in acid–base electrodialysis. Journal of Membrane Science, 182(1-2), pp.13-28. + +Wilhelm, F.G., Van Der Vegt, N.F.A., Strathmann, H. and Wessling, M., 2002. Comparison of bipolar membranes by means of chronopotentiometry. Journal of membrane science, 199(1-2), pp.177-190. \ No newline at end of file diff --git a/docs/technical_reference/unit_models/index.rst b/docs/technical_reference/unit_models/index.rst index 2015682160..3669c6d20f 100644 --- a/docs/technical_reference/unit_models/index.rst +++ b/docs/technical_reference/unit_models/index.rst @@ -6,6 +6,7 @@ Unit Models anaerobic_digester aeration_tank + bipolar_electrodialysis_0D boron_removal clarifier coag_floc_model diff --git a/watertap/costing/unit_models/bipolar_electrodialysis.py b/watertap/costing/unit_models/bipolar_electrodialysis.py new file mode 100644 index 0000000000..e70e83b9bf --- /dev/null +++ b/watertap/costing/unit_models/bipolar_electrodialysis.py @@ -0,0 +1,151 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# + +import pyomo.environ as pyo +from watertap.costing.util import ( + register_costing_parameter_block, + cost_rectifier, + make_capital_cost_var, + make_fixed_operating_cost_var, +) + + +def build_bipolar_electrodialysis_cost_param_block(blk): + # The following costing itemization and values are referenced to "Desalination 452 (2019) 265–278" + blk.membrane_capital_cost = pyo.Var( + initialize=160, + doc="Membrane and capital costs in [US$/m^2-membrane-area]", + units=pyo.units.USD_2018 / (pyo.units.meter**2), + ) + + blk.factor_membrane_replacement = pyo.Var( + initialize=0.2, + doc="Membrane and equipment (other stack components) housing replacement factor, equal to 1/lifetime.", + units=pyo.units.year**-1, + ) + + blk.stack_electrode_capital_cost = pyo.Var( + initialize=2100, + doc="Electrode cost in [US$/m^2-electrode-area] ", + units=pyo.units.USD_2018 / (pyo.units.meter**2), + ) + + blk.factor_stack_electrode_replacement = pyo.Var( + initialize=0.2, + doc="Stack and electrode replacement factor, equal to 1/lifetime.", + units=pyo.units.year**-1, + ) + + +@register_costing_parameter_block( + build_rule=build_bipolar_electrodialysis_cost_param_block, + parameter_block_name="bipolar_electrodialysis", +) +def cost_bipolar_electrodialysis( + blk, + cost_electricity_flow=True, + has_rectifier=False, +): + """ + Function for costing the bipolar electrodialysis unit + + Args: + cost_electricity_flow (:obj:`bool`, optional): Option for including the + costing of electricity. Defaults to True. + has_rectifier (:obj:`bool`, optional): Option for including a rectifier. + Defaults to False. + """ + t0 = blk.flowsheet().time.first() + + # Changed this to grab power from performance table which is identified + # by same key regardless of whether the Electrodialysis unit is 0D or 1D + if cost_electricity_flow: + if not has_rectifier: + blk.costing_package.cost_flow( + pyo.units.convert( + blk.unit_model.get_power_electrical(t0), + to_units=pyo.units.kW, + ), + "electricity", + ) + else: + power = blk.unit_model.get_power_electrical(blk.flowsheet().time.first()) + cost_rectifier(blk, power=power, ac_dc_conversion_efficiency=0.9) + + cost_bipolar_electrodialysis_stack(blk) + + +def cost_bipolar_electrodialysis_stack(blk): + """ + Generic function for costing the stack in an electrodialysis unit. + Assumes the unit_model has a `cell_num`, `cell_width`, and `cell_length` + set of variables used to size the total membrane area. + + """ + make_capital_cost_var(blk) + make_fixed_operating_cost_var(blk) + blk.costing_package.add_cost_factor(blk, "TIC") + if blk.find_component("capital_cost_rectifier") is not None: + blk.capital_cost_constraint = pyo.Constraint( + expr=blk.capital_cost + == blk.cost_factor + * ( + pyo.units.convert( + blk.costing_package.bipolar_electrodialysis.membrane_capital_cost + * ( + 2 + * blk.unit_model.cell_num + * blk.unit_model.cell_width + * blk.unit_model.cell_length + ) + + blk.costing_package.bipolar_electrodialysis.stack_electrode_capital_cost + * (2 * blk.unit_model.cell_width * blk.unit_model.cell_length), + to_units=blk.costing_package.base_currency, + ) + + blk.capital_cost_rectifier + ) + ) + else: + blk.capital_cost_constraint = pyo.Constraint( + expr=blk.capital_cost + == blk.cost_factor + * pyo.units.convert( + blk.costing_package.bipolar_electrodialysis.membrane_capital_cost + * ( + 2 + * blk.unit_model.cell_num + * blk.unit_model.cell_width + * blk.unit_model.cell_length + ) + + blk.costing_package.bipolar_electrodialysis.stack_electrode_capital_cost + * (2 * blk.unit_model.cell_width * blk.unit_model.cell_length), + to_units=blk.costing_package.base_currency, + ) + ) + blk.fixed_operating_cost_constraint = pyo.Constraint( + expr=blk.fixed_operating_cost + == pyo.units.convert( + blk.costing_package.bipolar_electrodialysis.factor_membrane_replacement + * blk.costing_package.bipolar_electrodialysis.membrane_capital_cost + * ( + 2 + * blk.unit_model.cell_num + * blk.unit_model.cell_width + * blk.unit_model.cell_length + ) + + blk.costing_package.bipolar_electrodialysis.factor_stack_electrode_replacement + * blk.costing_package.bipolar_electrodialysis.stack_electrode_capital_cost + * (2 * blk.unit_model.cell_width * blk.unit_model.cell_length), + to_units=blk.costing_package.base_currency + / blk.costing_package.base_period, + ) + ) diff --git a/watertap/unit_models/Bipolar_Electrodialysis_0D.py b/watertap/unit_models/Bipolar_Electrodialysis_0D.py new file mode 100644 index 0000000000..a0e5d41be6 --- /dev/null +++ b/watertap/unit_models/Bipolar_Electrodialysis_0D.py @@ -0,0 +1,2343 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# +import math + +# Import Pyomo libraries +from pyomo.environ import ( + Set, + Var, + check_optimal_termination, + Param, + Suffix, + NonNegativeReals, + value, + Constraint, + sqrt, + units as pyunits, +) +from pyomo.common.config import Bool, ConfigBlock, ConfigValue, In + +# Import IDAES cores +from idaes.core import ( + declare_process_block_class, + MaterialBalanceType, + EnergyBalanceType, + MomentumBalanceType, + UnitModelBlockData, + useDefault, +) +from idaes.core.util.misc import add_object_reference +from watertap.core.solvers import get_solver +from idaes.core.util.tables import create_stream_table_dataframe +from idaes.core.util.config import is_physical_parameter_block +from idaes.core.util.math import smooth_min + +from idaes.core.util.exceptions import ConfigurationError, InitializationError + +import idaes.core.util.scaling as iscale +import idaes.logger as idaeslog +from idaes.core.util.constants import Constants +from enum import Enum + +from watertap.core import ControlVolume0DBlock, InitializationMixin +from watertap.costing.unit_models.bipolar_electrodialysis import ( + cost_bipolar_electrodialysis, +) + +__author__ = " Johnson Dhanasekaran, Xiangyu Bi, Austin Ladshaw, Kejia Hu" + +_log = idaeslog.getLogger(__name__) + + +class LimitingCurrentDensitybpemMethod(Enum): + InitialValue = 0 + Empirical = 1 + + +class LimitingpotentialMethod(Enum): + InitialValue = 0 + Empirical = 1 + + +class ElectricalOperationMode(Enum): + Constant_Current = 0 + Constant_Voltage = 1 + + +class PressureDropMethod(Enum): + none = 0 + experimental = 1 + Darcy_Weisbach = 2 + + +class FrictionFactorMethod(Enum): + fixed = 0 + Gurreri = 1 + Kuroda = 2 + + +class HydraulicDiameterMethod(Enum): + fixed = 0 + spacer_specific_area_known = 1 + conventional = 2 + + +# Name of the unit model +@declare_process_block_class("Bipolar_Electrodialysis_0D") +class BipolarElectrodialysis0DData(InitializationMixin, UnitModelBlockData): + """ + 0D Bipolar and Electrodialysis Model + """ + + # CONFIG are options for the unit model + CONFIG = ConfigBlock() # + + CONFIG.declare( + "dynamic", + ConfigValue( + domain=In([False]), + default=False, + description="Dynamic model flag - must be False", + doc="""Indicates whether this model will be dynamic or not, + **default** = False. The filtration unit does not support dynamic + behavior, thus this must be False.""", + ), + ) + + CONFIG.declare( + "has_holdup", + ConfigValue( + default=False, + domain=In([False]), + description="Holdup construction flag - must be False", + doc="""Indicates whether holdup terms should be constructed or not. + **default** - False. The filtration unit does not have defined volume, thus + this must be False.""", + ), + ) + CONFIG.declare( + "has_pressure_change", + ConfigValue( + default=False, + domain=In([True, False]), + description="Pressure change term construction flag", + doc="""Indicates whether terms for pressure change should be + constructed, + **default** - False. + **Valid values:** { + **True** - include pressure change terms, + **False** - exclude pressure change terms.}""", + ), + ) + CONFIG.declare( + "pressure_drop_method", + ConfigValue( + default=PressureDropMethod.none, + domain=In(PressureDropMethod), + description="Method to calculate the frictional pressure drop in electrodialysis channels", + doc=""" + **default** - ``PressureDropMethod.none`` + + .. csv-table:: + :header: "Configuration Options", "Description" + + "``PressureDropMethod.none``", "The frictional pressure drop is neglected." + "``PressureDropMethod.experimental``", "The pressure drop is calculated by an experimental data as pressure drop per unit lenght." + "``PressureDropMethod.Darcy_Weisbach``", "The pressure drop is calculated by the Darcy-Weisbach equation." + """, + ), + ) + CONFIG.declare( + "friction_factor_method", + ConfigValue( + default=FrictionFactorMethod.fixed, + domain=In(FrictionFactorMethod), + description="Method to calculate the Darcy's friction factor", + doc=""" + **default** - ``FrictionFactorMethod.fixed`` + + .. csv-table:: + :header: "Configuration Options", "Description" + + "``FrictionFactorMethod.fixed``", "Friction factor is fixed by users" + "``FrictionFactorMethod.Gurreri``", "Friction factor evaluated based on Gurreri's work" + "``FrictionFactorMethod.Kuroda``", "Friction factor evaluated based on Kuroda's work" + """, + ), + ) + + CONFIG.declare( + "hydraulic_diameter_method", + ConfigValue( + default=HydraulicDiameterMethod.conventional, + domain=In(HydraulicDiameterMethod), + description="Method to calculate the hydraulic diameter for a rectangular channel in ED", + doc=""" + **default** - ``HydraulicDiameterMethod.conventional`` + + .. csv-table:: + :header: "Configuration Options", "Description" + + "``HydraulicDiameterMethod.fixed``", "Hydraulic diameter is fixed by users" + "``HydraulicDiameterMethod.conventional``", "Conventional method for a rectangular channel with spacer porosity considered" + "``HydraulicDiameterMethod.spacer_specific_area_known``", "A method for spacer-filled channel requiring the spacer specific area data" + """, + ), + ) + + CONFIG.declare( + "limiting_current_density_method_bpem", + ConfigValue( + default=LimitingCurrentDensitybpemMethod.InitialValue, + domain=In(LimitingCurrentDensitybpemMethod), + description="Configuration for method to compute the limiting current density across the bipolar membrane", + doc=""" + **default** - ``LimitingCurrentDensitybpemMethod.InitialValue`` + + .. csv-table:: + :header: "Configuration Options", "Description" + + "``LimitingCurrentDensitybpemMethod.InitialValue``", "Limiting current is calculated from a single initial value given by the user." + "``LimitingCurrentDensitybpemMethod.Empirical``", "Limiting current density is calculated from the empirical relationship" + """, + ), + ) + + CONFIG.declare( + "has_catalyst", + ConfigValue( + default=False, + domain=Bool, + description="""Catalyst action on water spliting, + **default** - False.""", + ), + ) + + CONFIG.declare( + "limiting_potential_method_bpem", + ConfigValue( + default=LimitingpotentialMethod.InitialValue, + domain=In(LimitingpotentialMethod), + description="Configuration for method to compute the limiting potential in bipolar membrane", + doc=""" + **default** - ``LimitingpotentialMethod.InitialValue`` + + .. csv-table:: + :header: "Configuration Options", "Description" + + "``LimitingpotentialMethod.InitialValue``", "Limiting current is calculated from a initial value given by the user." + "``LimitingpotentialMethod.Empirical``", "Limiting current density is caculated from the empirical equation" + """, + ), + ) + + CONFIG.declare( + "limiting_current_density_bpem_data", + ConfigValue( + default=0.5, + description="Limiting current density data input for bipolar membrane", + ), + ) + + CONFIG.declare( + "limiting_potential_data", + ConfigValue( + default=0.5, + description="Limiting potential of the bipolar membrane input", + ), + ) + + CONFIG.declare( + "material_balance_type", + ConfigValue( + default=MaterialBalanceType.useDefault, + domain=In(MaterialBalanceType), + description="Material balance construction flag", + doc="""Indicates what type of mass balance should be constructed, + **default** - MaterialBalanceType.useDefault. + **Valid values:** { + **MaterialBalanceType.useDefault - refer to property package for default + balance type + **MaterialBalanceType.none** - exclude material balances, + **MaterialBalanceType.componentPhase** - use phase component balances, + **MaterialBalanceType.componentTotal** - use total component balances, + **MaterialBalanceType.elementTotal** - use total element balances, + **MaterialBalanceType.total** - use total material balance.}""", + ), + ) + + CONFIG.declare( + "is_isothermal", + ConfigValue( + default=True, + domain=Bool, + description="""Assume isothermal conditions for control volume(s); energy_balance_type must be EnergyBalanceType.none, + **default** - True.""", + ), + ) + + CONFIG.declare( + "energy_balance_type", + ConfigValue( + default=EnergyBalanceType.none, + domain=In(EnergyBalanceType), + description="Energy balance construction flag", + doc="""Indicates what type of energy balance should be constructed, + **default** - EnergyBalanceType.none. + **Valid values:** { + **EnergyBalanceType.useDefault - refer to property package for default + balance type + **EnergyBalanceType.none** - exclude energy balances, + **EnergyBalanceType.enthalpyTotal** - single enthalpy balance for material, + **EnergyBalanceType.enthalpyPhase** - enthalpy balances for each phase, + **EnergyBalanceType.energyTotal** - single energy balance for material, + **EnergyBalanceType.energyPhase** - energy balances for each phase.}""", + ), + ) + + CONFIG.declare( + "momentum_balance_type", + ConfigValue( + default=MomentumBalanceType.pressureTotal, + domain=In(MomentumBalanceType), + description="Momentum balance construction flag", + doc="""Indicates what type of momentum balance should be constructed, + **default** - MomentumBalanceType.pressureTotal. + **Valid values:** { + **MomentumBalanceType.none** - exclude momentum balances, + **MomentumBalanceType.pressureTotal** - single pressure balance for material, + **MomentumBalanceType.pressurePhase** - pressure balances for each phase, + **MomentumBalanceType.momentumTotal** - single momentum balance for material, + **MomentumBalanceType.momentumPhase** - momentum balances for each phase.}""", + ), + ) + + CONFIG.declare( + "property_package", + ConfigValue( + default=useDefault, + domain=is_physical_parameter_block, + description="Property package to use for control volume", + doc="""Property parameter object used to define property calculations, + **default** - useDefault. + **Valid values:** { + **useDefault** - use default package from parent model or flowsheet, + **PhysicalParameterObject** - a PhysicalParameterBlock object.}""", + ), + ) + + CONFIG.declare( + "property_package_args", + ConfigBlock( + implicit=True, + description="Arguments to use for constructing property packages", + doc="""A ConfigBlock with arguments to be passed to a property block(s) + and used when constructing these, + **default** - None. + **Valid values:** { + see property package for documentation.}""", + ), + ) + + def _validate_config(self): + if ( + self.config.is_isothermal + and self.config.energy_balance_type != EnergyBalanceType.none + ): + raise ConfigurationError( + "If the isothermal assumption is used then the energy balance type must be none" + ) + + def build(self): + # build always starts by calling super().build() + # This triggers a lot of boilerplate in the background for you + super().build() + # this creates blank scaling factors, which are populated later + self.scaling_factor = Suffix(direction=Suffix.EXPORT) + + # Check configs for errors + self._validate_config() + + # Create essential sets. + self.membrane_set = Set(initialize=["bpem"]) + + add_object_reference(self, "ion_set", self.config.property_package.ion_set) + + add_object_reference( + self, "cation_set", self.config.property_package.cation_set + ) + add_object_reference(self, "anion_set", self.config.property_package.anion_set) + # Create unit model parameters and vars + self.water_density = Param( + initialize=1000, + units=pyunits.kg * pyunits.m**-3, + doc="density of water", + ) + + self.cell_num = Var( + initialize=1, + domain=NonNegativeReals, + bounds=(1, 10000), + units=pyunits.dimensionless, + doc="cell pair number in a stack", + ) + + # electrodialysis cell dimensional properties + self.cell_width = Var( + initialize=0.1, + bounds=(1e-3, 1e3), + units=pyunits.meter, + doc="The width of the electrodialysis cell, denoted as b in the model description", + ) + self.cell_length = Var( + initialize=0.5, + bounds=(1e-3, 1e2), + units=pyunits.meter, + doc="The length of the electrodialysis cell, denoted as l in the model description", + ) + self.channel_height = Var( + initialize=0.0001, + units=pyunits.meter, + doc="The distance between the consecutive aem and cem", + ) + self.spacer_porosity = Var( + initialize=0.7, + bounds=(0.01, 1), + units=pyunits.dimensionless, + doc='The prosity of spacer in the ED channels. This is also referred to elsewhere as "void fraction" or "volume parameters"', + ) + + # Material and Operational properties + self.membrane_thickness = Var( + self.membrane_set, + initialize=0.0001, + bounds=(1e-6, 1e-1), + units=pyunits.meter, + doc="Membrane thickness", + ) + self.solute_diffusivity_membrane = Var( + self.membrane_set, + self.ion_set | self.config.property_package.solute_set, + initialize=1e-10, + bounds=(0.0, 1e-6), + units=pyunits.meter**2 * pyunits.second**-1, + doc="Solute (ionic and neutral) diffusivity in the membrane phase", + ) + self.ion_trans_number_membrane = Var( + self.membrane_set, + self.ion_set, + bounds=(0, 1), + units=pyunits.dimensionless, + doc="Ion transference number in the membrane phase", + ) + self.water_trans_number_membrane = Var( + self.membrane_set, + initialize=5, + bounds=(0, 50), + units=pyunits.dimensionless, + doc="Transference number of water in membranes", + ) + self.water_permeability_membrane = Var( + self.membrane_set, + initialize=1e-14, + units=pyunits.meter * pyunits.second**-1 * pyunits.pascal**-1, + doc="Water permeability coefficient", + ) + self.membrane_areal_resistance = Var( + initialize=2e-4, + bounds=(0, 1), + units=pyunits.ohm * pyunits.meter**2, + doc="Surface resistance of membrane", + ) + self.electrodes_resistance = Var( + initialize=0, + bounds=(0, 100), + domain=NonNegativeReals, + units=pyunits.ohm * pyunits.meter**2, + doc="areal resistance of TWO electrode compartments of a stack", + ) + self.current = Var( + self.flowsheet().time, + initialize=1, + bounds=(0, 1e6), + units=pyunits.amp, + doc="Current across a cell-pair or stack", + ) + self.voltage = Var( + self.flowsheet().time, + initialize=100, + bounds=(0, 1000), + units=pyunits.volt, + doc="Voltage across a stack, declared under the 'Constant Voltage' mode only", + ) + self.current_utilization = Var( + initialize=1, + bounds=(0, 1), + units=pyunits.dimensionless, + doc="The current utilization including water electro-osmosis and ion diffusion", + ) + self.shadow_factor = Var( + initialize=1, + bounds=(0, 1), + units=pyunits.dimensionless, + doc="The reduction in area due to limited cross-section available for flow", + ) + + # Performance metrics + self.power_electrical = Var( + self.flowsheet().time, + initialize=1, + bounds=(0, 12100), + domain=NonNegativeReals, + units=pyunits.watt, + doc="Electrical power consumption of a stack", + ) + self.specific_power_electrical = Var( + self.flowsheet().time, + initialize=10, + bounds=(0, 1000), + domain=NonNegativeReals, + units=pyunits.kW * pyunits.hour * pyunits.meter**-3, + doc="acidic-volume-flow-rate-specific electrical power consumption", + ) + self.acid_produced = Var( + initialize=55 * 1e3, + bounds=(0, 1e6), + units=pyunits.kg * pyunits.second**-1, + doc="Acid prodcued", + ) + self.base_produced = Var( + initialize=55 * 1e3, + bounds=(0, 1e6), + units=pyunits.kg * pyunits.second**-1, + doc="Base prodcued", + ) + self.velocity_basic = Var( + self.flowsheet().time, + initialize=0.01, + units=pyunits.meter * pyunits.second**-1, + doc="Linear velocity of flow in the base channel of the bipolar membrane", + ) + self.velocity_acidic = Var( + self.flowsheet().time, + initialize=0.01, + units=pyunits.meter * pyunits.second**-1, + doc="Linear velocity of flow in the acid channel of the bipolar membrane", + ) + # Parameters for bipolar membrane operation + self.elec_field_non_dim = Var( + self.flowsheet().time, + initialize=1, + # bounds=(0, 1e32), + units=pyunits.dimensionless, + doc="Limiting current density across the bipolar membrane as a function of the normalized length", + ) + self.relative_permittivity = Var( + self.membrane_set, + initialize=30, + bounds=(1, 80), + domain=NonNegativeReals, + units=pyunits.dimensionless, + doc="Relative permittivity", + ) + self.membrane_fixed_charge = Var( + self.membrane_set, + initialize=1.5e3, + bounds=(1e-1, 1e5), + units=pyunits.mole * pyunits.meter**-3, + doc="Membrane fixed charge", + ) + self.kr = Var( + self.membrane_set, + initialize=1.33 * 10**11, + bounds=(1e-6, 1e16), + units=pyunits.L * pyunits.mole**-1 * pyunits.second**-1, + doc="Re-association rate constant", + ) + self.k2_zero = Var( + self.membrane_set, + initialize=2 * 10**-5, + bounds=(1e-10, 1e2), + units=pyunits.second**-1, + doc="Dissociation rate constant at no electric field", + ) + self.salt_conc_aem = Var( + self.membrane_set, + initialize=1e3, + bounds=(1e-8, 1e6), + units=pyunits.mole * pyunits.meter**-3, + doc="Salt concentration on the base channel of the bipolar membrane", + ) + self.salt_conc_cem = Var( + self.membrane_set, + initialize=1e3, + bounds=(1e-6, 1e4), + units=pyunits.mole * pyunits.meter**-3, + doc="Salt concentration on the acid channel of the bipolar membrane", + ) + self.diffus_mass = Var( + # self.membrane_set, + initialize=2e-9, + bounds=(1e-16, 1e-6), + units=pyunits.meter**2 * pyunits.second**-1, + doc="The mass diffusivity of the solute as molecules (not individual ions)", + ) + self.conc_water = Var( + self.membrane_set, + initialize=55 * 1e3, + bounds=(1e-2, 1e6), + units=pyunits.mole * pyunits.meter**-3, + doc="Concentration of water within the channel", + ) + + # Limiting operation quantity + self.current_dens_lim_bpem = Var( + self.flowsheet().time, + initialize=1e2, + bounds=(0, 1e5), + units=pyunits.amp * pyunits.meter**-2, + doc="Limiting current density across the bipolar membrane", + ) + + # Fluxes Vars for constructing mass transfer terms + self.generation_cem_flux_in = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_in of a component generated by water splitting on the acid channel of the bipolar membrane", + ) + self.generation_cem_flux_out = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_in of a component generated by water splitting on the acid channel of the bipolar membrane", + ) + self.generation_aem_flux_in = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_in of a component generated by water splitting on the base channel of the bipolar membrane", + ) + self.generation_aem_flux_out = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_in of a component generated by water splitting on the base channel of the bipolar membrane", + ) + self.elec_migration_bpem_flux_in = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_in of a component across the membrane driven by electrical migration across the bipolar membrane", + ) + self.elec_migration_bpem_flux_out = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_out of a component across the membrane driven by electrical migration across the bipolar membrane", + ) + self.nonelec_bpem_flux_in = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_in of a component across the membrane driven by non-electrical forces across the bipolar membrane", + ) + self.nonelec_bpem_flux_out = Var( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Molar flux_out of a component across the membrane driven by non-electrical forces across the bipolar membrane", + ) + + # Build control volume for the base channel of the bipolar channel + self.basic = ControlVolume0DBlock( + dynamic=False, + has_holdup=False, + property_package=self.config.property_package, + property_package_args=self.config.property_package_args, + ) + self.basic.add_state_blocks(has_phase_equilibrium=False) + self.basic.add_material_balances( + balance_type=self.config.material_balance_type, has_mass_transfer=True + ) + self.basic.add_energy_balances( + balance_type=self.config.energy_balance_type, + has_enthalpy_transfer=False, + ) + + if self.config.is_isothermal: + self.basic.add_isothermal_assumption() + self.basic.add_momentum_balances( + balance_type=self.config.momentum_balance_type, + has_pressure_change=self.config.has_pressure_change, + ) + # Build control volume for the acid channel of the bipolar membrane channel + self.acidic = ControlVolume0DBlock( + dynamic=False, + has_holdup=False, + property_package=self.config.property_package, + property_package_args=self.config.property_package_args, + ) + self.acidic.add_state_blocks(has_phase_equilibrium=False) + self.acidic.add_material_balances( + balance_type=self.config.material_balance_type, has_mass_transfer=True + ) + self.acidic.add_energy_balances( + balance_type=self.config.energy_balance_type, + has_enthalpy_transfer=False, + ) + + if self.config.is_isothermal: + self.acidic.add_isothermal_assumption() + self.acidic.add_momentum_balances( + balance_type=self.config.momentum_balance_type, + has_pressure_change=self.config.has_pressure_change, + ) + + # den_mass and visc_d in acidic and basic channels are the same + add_object_reference( + self, "dens_mass", self.acidic.properties_in[0].dens_mass_phase["Liq"] + ) + add_object_reference( + self, "visc_d", self.acidic.properties_in[0].visc_d_phase["Liq"] + ) + + # Add ports (creates inlets and outlets for each channel) + self.add_inlet_port(name="inlet_basic", block=self.basic) + self.add_outlet_port(name="outlet_basic", block=self.basic) + self.add_inlet_port(name="inlet_acidic", block=self.acidic) + self.add_outlet_port(name="outlet_acidic", block=self.acidic) + + # extension options + if self.config.has_catalyst == True: + self._make_catalyst() + + if ( + not self.config.pressure_drop_method == PressureDropMethod.none + ) and self.config.has_pressure_change: + self._pressure_drop_calculation() + + @self.Constraint( + self.flowsheet().time, + doc="Pressure drop expression as calculated by the pressure drop data, " + "base channel of the bipolar membrane.", + ) + def eq_deltaP_basic(self, t): + return self.basic.deltaP[t] == -self.pressure_drop_total[t] + + @self.Constraint( + self.flowsheet().time, + doc="Pressure drop expression as calculated by the pressure drop data," + " acid channel of the bipolar membrane.", + ) + def eq_deltaP_acidic(self, t): + return self.acidic.deltaP[t] == -self.pressure_drop_total[t] + + elif self.config.pressure_drop_method == PressureDropMethod.none and ( + not self.config.has_pressure_change + ): + pass + else: + raise ConfigurationError( + "A valid (not none) pressure_drop_method and has_pressure_change being True " + "must be both used or unused at the same time. " + ) + + # Build Constraints + + @self.Constraint( + self.flowsheet().time, + doc="Calculate flow velocity in a single base (A.E.M side) channel of the bipolar membrane channel," + " based on the average of inlet and outlet", + ) + def eq_get_velocity_basic(self, t): + return self.velocity_basic[ + t + ] * self.cell_width * self.shadow_factor * self.channel_height * self.spacer_porosity * self.cell_num == 0.5 * ( + self.basic.properties_in[0].flow_vol_phase["Liq"] + + self.basic.properties_out[0].flow_vol_phase["Liq"] + ) + + @self.Constraint( + self.flowsheet().time, + doc="Calculate flow velocity in a single acid (C.E.M side) channel of the bipolar membrane channel," + " based on the average of inlet and outlet", + ) + def eq_get_velocity_acidic(self, t): + return self.velocity_acidic[ + t + ] * self.cell_width * self.shadow_factor * self.channel_height * self.spacer_porosity * self.cell_num == 0.5 * ( + self.acidic.properties_in[0].flow_vol_phase["Liq"] + + self.acidic.properties_out[0].flow_vol_phase["Liq"] + ) + + @self.Constraint( + self.flowsheet().time, + doc="Calculate limiting current density across the bipolar membrane", + ) + def eq_current_dens_lim_bpem(self, t): + if ( + self.config.limiting_current_density_method_bpem + == LimitingCurrentDensitybpemMethod.InitialValue + ): + return self.current_dens_lim_bpem[t] == ( + self.config.limiting_current_density_bpem_data + * pyunits.amp + * pyunits.meter**-2 + ) + elif ( + self.config.limiting_current_density_method_bpem + == LimitingCurrentDensitybpemMethod.Empirical + ): + return self.current_dens_lim_bpem[ + t + ] == self.diffus_mass * Constants.faraday_constant * ( + (self.salt_conc_aem["bpem"] + self.salt_conc_cem["bpem"]) * 0.5 + ) ** 2 / ( + self.membrane_thickness["bpem"] * self.membrane_fixed_charge["bpem"] + ) + + @self.Constraint( + self.flowsheet().time, + doc="Calculate the potential drops across the bipolar membrane", + ) + def eq_potential_barrier_bpem(self, t): + if self.config.has_catalyst: + return Constraint.Skip + else: + self.potential_barrier_bpem = Var( + self.flowsheet().time, + initialize=1, + bounds=(0, 5000), + units=pyunits.volt, + doc="Potential barrier across the depletion layer for water splitting to begin", + ) + + if ( + self.config.limiting_potential_method_bpem + == LimitingpotentialMethod.InitialValue + ): + return self.potential_barrier_bpem[t] == ( + self.config.limiting_potential_data * pyunits.volt + ) + + elif ( + self.config.limiting_potential_method_bpem + == LimitingpotentialMethod.Empirical + ): + # [H+][OH-] concentration + kw = 10**-8 * pyunits.mol**2 * pyunits.meter**-6 + + # Fraction of threshold of limiting current: currently 0.1 i_lim + frac = 1 * 10**-1 + # Dimensional pre-factor to evaulate non-dimensional electric field + const = 0.0936 * pyunits.K**2 * pyunits.volt**-1 * pyunits.meter + + @self.Constraint( + self.flowsheet().time, + doc="Calculate the non-dimensional potential drop", + ) + def eq_potential_barrier_bpem_non_dim(self, t): + # [y2, qty_une, qty_deux, qty_trois] = dat + terms = 40 + matrx = 0 + for indx in range(terms): + # rev_indx = terms - indx - 1 + matrx += ( + 2**indx + * self.elec_field_non_dim[t] ** indx + / (math.factorial(indx) * math.factorial(indx + 1)) + ) + + matrx *= self.k2_zero["bpem"] * self.conc_water["bpem"] + matrx += ( + -pyunits.convert( + self.kr["bpem"], + to_units=pyunits.meter**3 + * pyunits.mole**-1 + * pyunits.second**-1, + ) + * kw + ) + return ( + Constants.vacuum_electric_permittivity + * self.relative_permittivity["bpem"] ** 2 + * self.basic.properties_in[t].temperature ** 2 + * Constants.avogadro_number + * Constants.elemental_charge + ) / ( + const + * Constants.faraday_constant + * self.membrane_fixed_charge["bpem"] + ) * matrx * self.elec_field_non_dim[ + t + ] == self.current_dens_lim_bpem[ + t + ] * frac + + # Dimensional electric field + field_generated = ( + self.elec_field_non_dim[t] + * self.relative_permittivity["bpem"] + * self.basic.properties_in[t].temperature ** 2 + / const + ) + + # Depletion length at the junction of the bipolar membrane + lambda_depletion = ( + field_generated + * Constants.vacuum_electric_permittivity + * self.relative_permittivity["bpem"] + / ( + Constants.faraday_constant + * self.membrane_fixed_charge["bpem"] + ) + ) + + return ( + self.potential_barrier_bpem[t] + == field_generated * lambda_depletion + ) + + else: + self.potential_barrier_bpem[t].fix(0 * pyunits.volt) + return Constraint.Skip + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + doc="Current-Voltage relationship", + ) + def eq_current_voltage_relation(self, t, p): + + total_areal_resistance = ( + self.membrane_areal_resistance + + self.channel_height + * ( + 0.5**-1 + * ( + self.basic.properties_in[t].elec_cond_phase["Liq"] + + self.basic.properties_out[t].elec_cond_phase["Liq"] + ) + ** -1 + + 0.5**-1 + * ( + self.acidic.properties_in[t].elec_cond_phase["Liq"] + + self.acidic.properties_out[t].elec_cond_phase["Liq"] + ) + ** -1 + ) + ) * self.cell_num + self.electrodes_resistance + # the average conductivity of each channel's inlet and outlet is taken to represent that of the entire channel + + if self.config.has_catalyst: + voltage_membrane_drop = self.potential_membrane_bpem[t] + + @self.Constraint( + self.flowsheet().time, + doc="Calculate total current generated via catalyst action", + ) + def eq_current_relationship(self, t): + return self.current[t] == ( + self.current_dens_lim_bpem[t] + + self.flux_splitting[t] * Constants.faraday_constant + ) * (self.cell_width * self.shadow_factor * self.cell_length) + + else: + voltage_membrane_drop = self.potential_barrier_bpem[t] + + return ( + self.current[t] + * (self.cell_width * self.shadow_factor * self.cell_length) ** -1 + * total_areal_resistance + + voltage_membrane_drop * self.cell_num + == self.voltage[t] + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for water splitting acid channel (C.E.M Side) of bipolar membrane flux_in", + ) + def eq_generation_cem_flux_in(self, t, p, j): + if j == "H_+": + if self.config.has_catalyst == True: + return ( + self.generation_cem_flux_in[t, p, j] == self.flux_splitting[t] + ) + + else: + return self.generation_cem_flux_in[t, p, j] == ( + -smooth_min( + -( + self.current[t] / pyunits.amp + - self.current_dens_lim_bpem[t] + * self.cell_width + * self.shadow_factor + * self.cell_length + / pyunits.amp + ), + 0, + ) + * pyunits.amp + / (self.cell_width * self.shadow_factor * self.cell_length) + ) / (Constants.faraday_constant) + + else: + if j == "H2O": + if self.config.has_catalyst == True: + return ( + self.generation_cem_flux_in[t, p, j] + == -0.5 * self.flux_splitting[t] + ) + else: + return self.generation_cem_flux_in[t, p, j] == ( + 0.5 + * smooth_min( + -( + self.current[t] / pyunits.amp + - self.current_dens_lim_bpem[t] + * self.cell_width + * self.shadow_factor + * self.cell_length + / pyunits.amp + ), + 0, + ) + * pyunits.amp + / (self.cell_width * self.shadow_factor * self.cell_length) + ) / (Constants.faraday_constant) + + else: + return ( + self.generation_cem_flux_in[t, p, j] + == 0 * pyunits.mol * pyunits.meter**-2 * pyunits.s**-1 + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for water splitting base channel (A.E.M Side) of bipolar membrane flux_in", + ) + def eq_generation_aem_flux_in(self, t, p, j): + if j == "OH_-": + if self.config.has_catalyst == True: + return ( + self.generation_aem_flux_in[t, p, j] == self.flux_splitting[t] + ) + + else: + return self.generation_aem_flux_in[t, p, j] == ( + -smooth_min( + -( + self.current[t] / pyunits.amp + - self.current_dens_lim_bpem[t] + * self.cell_width + * self.shadow_factor + * self.cell_length + / pyunits.amp + ), + 0, + ) + * pyunits.amp + / (self.cell_width * self.shadow_factor * self.cell_length) + ) / (Constants.faraday_constant) + + else: + if j == "H2O": + if self.config.has_catalyst == True: + return ( + self.generation_aem_flux_in[t, p, j] + == -0.5 * self.flux_splitting[t] + ) + + else: + return self.generation_aem_flux_in[t, p, j] == ( + 0.5 + * smooth_min( + -( + self.current[t] / pyunits.amp + - self.current_dens_lim_bpem[t] + * self.cell_width + * self.shadow_factor + * self.cell_length + / pyunits.amp + ), + 0, + ) + * pyunits.amp + / (self.cell_width * self.shadow_factor * self.cell_length) + ) / (Constants.faraday_constant) + + else: + return ( + self.generation_aem_flux_in[t, p, j] + == 0 * pyunits.mol * pyunits.meter**-2 * pyunits.s**-1 + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for water splitting acid channel (C.E.M Side) of bipolar membrane flux_out", + ) + def eq_generation_cem_flux_out(self, t, p, j): + return ( + self.generation_cem_flux_in[t, p, j] + == self.generation_cem_flux_out[t, p, j] + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for water splitting base channel (A.E.M Side) of bipolar membrane flux_out", + ) + def eq_generation_aem_flux_out(self, t, p, j): + return ( + self.generation_aem_flux_in[t, p, j] + == self.generation_aem_flux_out[t, p, j] + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for electrical migration across the bipolar membrane flux_in", + ) + def eq_elec_migration_bpem_flux_in(self, t, p, j): + if j == "H2O": + return self.elec_migration_bpem_flux_in[t, p, j] == ( + self.water_trans_number_membrane["bpem"] + ) * ( + self.current[t] + / (self.cell_width * self.shadow_factor * self.cell_length) + / Constants.faraday_constant + ) + + elif j in self.ion_set: + if not (j == "H_+" or j == "OH_-"): + if self.config.has_catalyst == False: + + return self.elec_migration_bpem_flux_in[t, p, j] == ( + self.ion_trans_number_membrane["bpem", j] + ) * ( + self.current_utilization + * smooth_min( + self.current[t] / pyunits.amp, + self.current_dens_lim_bpem[t] + * self.cell_width + * self.shadow_factor + * self.cell_length + / pyunits.amp, + ) + * pyunits.amp + / (self.cell_width * self.shadow_factor * self.cell_length) + ) / ( + self.config.property_package.charge_comp[j] + * Constants.faraday_constant + ) + else: + return self.elec_migration_bpem_flux_in[t, p, j] == ( + self.ion_trans_number_membrane["bpem", j] + ) * ( + self.current_utilization * self.current_dens_lim_bpem[t] + ) / ( + self.config.property_package.charge_comp[j] + * Constants.faraday_constant + ) + + else: + + self.elec_migration_bpem_flux_in[t, p, j].fix( + 0 * pyunits.mol * pyunits.m**-2 * pyunits.s**-1 + ) + return Constraint.Skip + else: + self.elec_migration_bpem_flux_in[t, p, j].fix( + 0 * pyunits.mol * pyunits.m**-2 * pyunits.s**-1 + ) + return Constraint.Skip + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for electrical migration across the bipolar membrane flux_out", + ) + def eq_elec_migration_bpem_flux_out(self, t, p, j): + return ( + self.elec_migration_bpem_flux_out[t, p, j] + == self.elec_migration_bpem_flux_in[t, p, j] + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for non-electrical flux across the bipolar membrane flux_in", + ) + def eq_nonelec_bpem_flux_in(self, t, p, j): + if j == "H2O": + return self.nonelec_bpem_flux_in[ + t, p, j + ] == self.water_density / self.config.property_package.mw_comp[j] * ( + self.water_permeability_membrane["bpem"] + ) * ( + self.basic.properties_in[t].pressure_osm_phase[p] + - self.acidic.properties_in[t].pressure_osm_phase[p] + ) + + else: + return ( + self.nonelec_bpem_flux_in[t, p, j] + == 0 * pyunits.mol * pyunits.m**-2 * pyunits.s**-1 + ) + + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Equation for non-electrical flux across the bipolar membrane flux_out", + ) + def eq_nonelec_bpem_flux_out(self, t, p, j): + if j == "H2O": + return self.nonelec_bpem_flux_out[ + t, p, j + ] == self.water_density / self.config.property_package.mw_comp[j] * ( + self.water_permeability_membrane["bpem"] + ) * ( + self.basic.properties_out[t].pressure_osm_phase[p] + - self.acidic.properties_out[t].pressure_osm_phase[p] + ) + + else: + return ( + self.nonelec_bpem_flux_out[t, p, j] + == 0 * pyunits.mol * pyunits.m**-2 * pyunits.s**-1 + ) + + # Add constraints for mass transfer terms (base channel/A.E.M Side of the bipolar membrane) + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Mass transfer term for the base channel (A.E.M Side) of the bipolar membrane", + ) + def eq_mass_transfer_term_basic(self, t, p, j): + return ( + self.basic.mass_transfer_term[t, p, j] + == 0.5 + * ( + self.generation_aem_flux_in[t, p, j] + + self.generation_aem_flux_out[t, p, j] + + self.elec_migration_bpem_flux_in[t, p, j] + + self.elec_migration_bpem_flux_out[t, p, j] + + self.nonelec_bpem_flux_in[t, p, j] + + self.nonelec_bpem_flux_out[t, p, j] + ) + * (self.cell_width * self.shadow_factor * self.cell_length) + * self.cell_num + ) + + # Add constraints for mass transfer terms (acid channel/C.E.M Side of the bipolar membrane) + @self.Constraint( + self.flowsheet().time, + self.config.property_package.phase_list, + self.config.property_package.component_list, + doc="Mass transfer term for the acid channel (C.E.M Side) of the bipolar membrane channel", + ) + def eq_mass_transfer_term_acidic(self, t, p, j): + return ( + self.acidic.mass_transfer_term[t, p, j] + == 0.5 + * ( + self.generation_cem_flux_in[t, p, j] + + self.generation_cem_flux_out[t, p, j] + + self.elec_migration_bpem_flux_in[t, p, j] + + self.elec_migration_bpem_flux_out[t, p, j] + + self.nonelec_bpem_flux_in[t, p, j] + + self.nonelec_bpem_flux_out[t, p, j] + ) + * (self.cell_width * self.shadow_factor * self.cell_length) + * self.cell_num + ) + + # Performance: acid/base produced + @self.Constraint( + self.flowsheet().time, + doc="Evaluate Base produced", + ) + def eq_product_basic(self, t): + conc_unit = 1 * pyunits.mole * pyunits.second**-1 + product_net_loc = 0 * pyunits.kg * pyunits.second**-1 + + for j in self.config.property_package.cation_set: + if not j == "H_+": + product_in_loc = smooth_min( + self.basic.properties_in[t].flow_mol_phase_comp["Liq", j] + / conc_unit, + self.basic.properties_in[t].flow_mol_phase_comp["Liq", "OH_-"] + / conc_unit + / self.config.property_package.charge_comp[j], + ) + + product_out_loc = smooth_min( + self.basic.properties_out[t].flow_mol_phase_comp["Liq", j] + / conc_unit, + self.basic.properties_out[t].flow_mol_phase_comp["Liq", "OH_-"] + / conc_unit + / self.config.property_package.charge_comp[j], + ) + + product_net_loc += ( + -1 + * smooth_min(product_in_loc - product_out_loc, 0) + * conc_unit + * ( + self.config.property_package.charge_comp[j] + * self.config.property_package.mw_comp["OH_-"] + + self.config.property_package.mw_comp[j] + ) + ) + + return self.base_produced == product_net_loc + + @self.Constraint( + self.flowsheet().time, + doc="Evaluate Acid produced", + ) + def eq_product_acidic(self, t): + conc_unit = 1 * pyunits.mole * pyunits.second**-1 + product_net_loc = 0 * pyunits.kg * pyunits.second**-1 + + for j in self.config.property_package.anion_set: + if not j == "OH_-": + product_in_loc = smooth_min( + self.acidic.properties_in[t].flow_mol_phase_comp["Liq", j] + / conc_unit, + self.acidic.properties_in[t].flow_mol_phase_comp["Liq", "H_+"] + / conc_unit + / (-self.config.property_package.charge_comp[j]), + ) + + product_out_loc = smooth_min( + self.acidic.properties_out[t].flow_mol_phase_comp["Liq", j] + / conc_unit, + self.acidic.properties_out[t].flow_mol_phase_comp["Liq", "H_+"] + / conc_unit + / (-self.config.property_package.charge_comp[j]), + ) + + product_net_loc += ( + -1 + * smooth_min(product_in_loc - product_out_loc, 0) + * conc_unit + * ( + (-self.config.property_package.charge_comp[j]) + * self.config.property_package.mw_comp["H_+"] + + self.config.property_package.mw_comp[j] + ) + ) + + return self.acid_produced == product_net_loc + + # Performance: Electrical + @self.Constraint( + self.flowsheet().time, + doc="Electrical power consumption of a stack", + ) + def eq_power_electrical(self, t): + return self.power_electrical[t] == self.current[t] * self.voltage[t] + + @self.Constraint( + self.flowsheet().time, + doc="Diluate_volume_flow_rate_specific electrical power consumption of a stack", + ) + def eq_specific_power_electrical(self, t): + return ( + pyunits.convert( + self.specific_power_electrical[t], + pyunits.watt * pyunits.second * pyunits.meter**-3, + ) + * self.acidic.properties_out[t].flow_vol_phase["Liq"] + == self.current[t] * self.voltage[t] + ) + + # Catalyst action: + def _make_catalyst(self): + + self.potential_membrane_bpem = Var( + self.flowsheet().time, + initialize=0.1, + bounds=(0, 1e8), + units=pyunits.volt, + doc="Potential drop across the depletion layer", + ) + self.flux_splitting = Var( + self.flowsheet().time, + initialize=1, + domain=NonNegativeReals, + # bounds=(0, 50000), + units=pyunits.mole * pyunits.meter**-2 * pyunits.second**-1, + doc="Flux generated", + ) + self.membrane_fixed_catalyst_aem = Var( + self.membrane_set, + initialize=5e3, + bounds=(1e-1, 1e5), + units=pyunits.mole * pyunits.meter**-3, + doc="Membrane fixed charge", + ) + self.membrane_fixed_catalyst_cem = Var( + self.membrane_set, + initialize=5e3, + bounds=(1e-1, 1e5), + units=pyunits.mole * pyunits.meter**-3, + doc="Membrane fixed charge", + ) + + self.k_a = Var( + initialize=1e-3, + bounds=(1e-15, 1e15), + units=pyunits.mole * pyunits.meter**-3, + doc="Membrane fixed charge", + ) + self.k_b = Var( + initialize=3e-2, + bounds=(1e-15, 1e15), + units=pyunits.mole * pyunits.meter**-3, + doc="Membrane fixed charge", + ) + + const = 0.0936 * pyunits.K**2 * pyunits.volt**-1 * pyunits.meter + + @self.Constraint( + self.flowsheet().time, + doc="Calculate the non-dimensional potential drop across the depletion region", + ) + def eq_potential_membrane_bpem_non_dim(self, t): + return self.elec_field_non_dim[t] == const * self.basic.properties_in[ + t + ].temperature ** -2 * self.relative_permittivity["bpem"] ** -1 * sqrt( + ( + Constants.faraday_constant + * self.membrane_fixed_charge["bpem"] + * self.potential_membrane_bpem[t] + ) + / ( + Constants.vacuum_electric_permittivity + * self.relative_permittivity["bpem"] + ) + ) + + @self.Constraint( + self.flowsheet().time, + doc="Calculate the potential barrier at limiting current across the bipolar membrane", + ) + def eq_flux_splitting(self, t): + terms = 40 + matrx = 0 + for indx in range(terms): + matrx += ( + 2**indx + * self.elec_field_non_dim[t] ** indx + / (math.factorial(indx) * math.factorial(indx + 1)) + ) + + matrx *= self.k2_zero["bpem"] * self.conc_water["bpem"] + matrx_a = matrx * self.membrane_fixed_catalyst_cem["bpem"] / self.k_a + matrx_b = matrx * self.membrane_fixed_catalyst_aem["bpem"] / self.k_b + return self.flux_splitting[t] == (matrx_a + matrx_b) * sqrt( + self.potential_membrane_bpem[t] + * Constants.vacuum_electric_permittivity + * self.relative_permittivity["bpem"] + / (Constants.faraday_constant * self.membrane_fixed_charge["bpem"]) + ) + + def _get_fluid_dimensionless_quantities(self): + # self.diffus_mass = Var( + # initialize=1e-9, + # bounds=(1e-16, 1e-6), + # units=pyunits.meter**2 * pyunits.second**-1, + # doc="The mass diffusivity of the solute as molecules (not individual ions)", + # ) + self.hydraulic_diameter = Var( + initialize=1e-3, + bounds=(0, None), + units=pyunits.meter, + doc="The hydraulic diameter of the channel", + ) + self.N_Re = Var( + initialize=50, + bounds=(0, None), + units=pyunits.dimensionless, + doc="Reynolds Number", + ) + # self.N_Sc = Var( + # initialize=2000, + # bounds=(0, None), + # units=pyunits.dimensionless, + # doc="Schmidt Number", + # ) + # self.N_Sh = Var( + # initialize=100, + # bounds=(0, None), + # units=pyunits.dimensionless, + # doc="Sherwood Number", + # ) + + if self.config.hydraulic_diameter_method == HydraulicDiameterMethod.fixed: + _log.warning("Do not forget to FIX the channel hydraulic diameter in [m]!") + else: + + @self.Constraint( + doc="To calculate hydraulic diameter", + ) + def eq_hydraulic_diameter(self): + if ( + self.config.hydraulic_diameter_method + == HydraulicDiameterMethod.conventional + ): + return ( + self.hydraulic_diameter + == 2 + * self.channel_height + * self.cell_width + * self.shadow_factor + * self.spacer_porosity + * (self.channel_height + self.cell_width * self.shadow_factor) + ** -1 + ) + else: + self.spacer_specific_area = Var( + initialize=1e4, + bounds=(0, None), + units=pyunits.meter**-1, + doc="The specific area of the channel", + ) + return ( + self.hydraulic_diameter + == 4 + * self.spacer_porosity + * ( + 2 * self.channel_height**-1 + + (1 - self.spacer_porosity) * self.spacer_specific_area + ) + ** -1 + ) + + @self.Constraint( + doc="To calculate Re", + ) + def eq_Re(self): + + return ( + self.N_Re + == self.dens_mass + * self.velocity_acidic[0] + * self.hydraulic_diameter + * self.visc_d**-1 + ) + + # @self.Constraint( + # doc="To calculate Sc", + # ) + # def eq_Sc(self): + # + # return self.N_Sc == self.visc_d * self.dens_mass**-1 * self.diffus_mass**-1 + # + # @self.Constraint( + # doc="To calculate Sh", + # ) + # def eq_Sh(self): + # + # return self.N_Sh == 0.29 * self.N_Re**0.5 * self.N_Sc**0.33 + + def _pressure_drop_calculation(self): + self.pressure_drop_total = Var( + self.flowsheet().time, + initialize=1e6, + units=pyunits.pascal, + doc="pressure drop over an entire ED stack", + ) + self.pressure_drop = Var( + self.flowsheet().time, + initialize=1e3, + units=pyunits.pascal * pyunits.meter**-1, + doc="pressure drop per unit of length", + ) + + if self.config.pressure_drop_method == PressureDropMethod.experimental: + _log.warning( + "Do not forget to FIX the experimental pressure drop value in [Pa/m]!" + ) + else: # PressureDropMethod.Darcy_Weisbach is used + # if not (self.config.has_Nernst_diffusion_layer == True): + self._get_fluid_dimensionless_quantities() + + self.friction_factor = Var( + initialize=10, + bounds=(0, None), + units=pyunits.dimensionless, + doc="friction factor of the channel fluid", + ) + + @self.Constraint( + self.flowsheet().time, + doc="To calculate pressure drop over an stack", + ) + def eq_pressure_drop_total(self, t): + return ( + self.pressure_drop_total[t] + == self.dens_mass + * self.friction_factor + * self.velocity_acidic[0] ** 2 + * 0.5 + * self.hydraulic_diameter**-1 + * self.cell_length + ) + + if self.config.friction_factor_method == FrictionFactorMethod.fixed: + _log.warning("Do not forget to FIX the Darcy's friction factor value!") + else: + + @self.Constraint( + doc="To calculate friction factor", + ) + def eq_friction_factor(self): + if ( + self.config.friction_factor_method + == FrictionFactorMethod.Gurreri + ): + return ( + self.friction_factor + == 4 * 50.6 * self.spacer_porosity**-7.06 * self.N_Re**-1 + ) + elif ( + self.config.friction_factor_method + == FrictionFactorMethod.Kuroda + ): + return ( + self.friction_factor + == 4 * 9.6 * self.spacer_porosity**-1 * self.N_Re**-0.5 + ) + + @self.Constraint( + self.flowsheet().time, + doc="To calculate total pressure drop over a stack", + ) + def eq_pressure_drop(self, t): + return ( + self.pressure_drop[t] + == self.pressure_drop_total[t] * self.cell_length**-1 + ) + + # initialize method + def initialize_build( + self, + state_args=None, + outlvl=idaeslog.NOTSET, + solver=None, + optarg=None, + ): + """ + General wrapper for pressure changer initialization routines + + Keyword Arguments: + state_args : a dict of arguments to be passed to the property + package(s) to provide an initial state for + initialization (see documentation of the specific + property package) (default = {}). + outlvl : sets output level of initialization routine + optarg : solver options dictionary object (default=None) + solver : str indicating which solver to use during + initialization (default = None) + + Returns: None + """ + init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit") + solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit") + # Set solver options + opt = get_solver(solver, optarg) + + # --------------------------------------------------------------------- + + # Set the outlet has the same intial condition of the inlet. + for k in self.keys(): + for j in self[k].config.property_package.component_list: + self[k].basic.properties_out[0].flow_mol_phase_comp["Liq", j] = value( + self[k].basic.properties_in[0].flow_mol_phase_comp["Liq", j] + ) + self[k].acidic.properties_out[0].flow_mol_phase_comp["Liq", j] = value( + self[k].acidic.properties_in[0].flow_mol_phase_comp["Liq", j] + ) + + # --------------------------------------------------------------------- + # Initialize concentrate_basic_side block + flags_basic = self.basic.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args, + hold_state=True, + ) + init_log.info_high("Initialization Step 2 Complete.") + # --------------------------------------------------------------------- + # Initialize concentrate_acidic_side block + flags_acidic = self.acidic.initialize( + outlvl=outlvl, + optarg=optarg, + solver=solver, + state_args=state_args, # inlet var + hold_state=True, + ) + init_log.info_high("Initialization Step 3 Complete.") + # --------------------------------------------------------------------- + # Solve unit + with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc: + res = opt.solve(self, tee=slc.tee) + init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res))) + # --------------------------------------------------------------------- + # Release state + self.basic.release_state(flags_basic, outlvl) + init_log.info("Initialization Complete: {}".format(idaeslog.condition(res))) + self.acidic.release_state(flags_acidic, outlvl) + init_log.info("Initialization Complete: {}".format(idaeslog.condition(res))) + + if not check_optimal_termination(res): + raise InitializationError(f"Unit model {self.name} failed to initialize") + + def calculate_scaling_factors(self): + super().calculate_scaling_factors() + # Scaling factors that user may setup + # The users are highly encouraged to provide scaling factors for assessable vars below. + # Not providing these vars will give a warning. + if ( + iscale.get_scaling_factor(self.solute_diffusivity_membrane, warning=True) + is None + ): + iscale.set_scaling_factor(self.solute_diffusivity_membrane, 1e10) + if iscale.get_scaling_factor(self.membrane_thickness, warning=True) is None: + iscale.set_scaling_factor(self.membrane_thickness, 1e4) + if ( + iscale.get_scaling_factor(self.water_permeability_membrane, warning=True) + is None + ): + iscale.set_scaling_factor(self.water_permeability_membrane, 1e14) + if iscale.get_scaling_factor(self.cell_num, warning=True) is None: + iscale.set_scaling_factor(self.cell_num, 0.1) + if iscale.get_scaling_factor(self.cell_length, warning=True) is None: + iscale.set_scaling_factor(self.cell_length, 1e1) + if iscale.get_scaling_factor(self.cell_width, warning=True) is None: + iscale.set_scaling_factor(self.cell_width, 1e2) + if iscale.get_scaling_factor(self.shadow_factor, warning=True) is None: + iscale.set_scaling_factor(self.shadow_factor, 1) + if iscale.get_scaling_factor(self.channel_height, warning=True) is None: + iscale.set_scaling_factor(self.channel_height, 1e5) + if iscale.get_scaling_factor(self.spacer_porosity, warning=True) is None: + iscale.set_scaling_factor(self.spacer_porosity, 1) + if ( + iscale.get_scaling_factor(self.membrane_areal_resistance, warning=True) + is None + ): + iscale.set_scaling_factor(self.membrane_areal_resistance, 1e5) + if iscale.get_scaling_factor(self.electrodes_resistance, warning=True) is None: + iscale.set_scaling_factor(self.electrodes_resistance, 1e1) + if iscale.get_scaling_factor(self.current, warning=True) is None: + iscale.set_scaling_factor(self.current, 1) + if iscale.get_scaling_factor(self.voltage, warning=True) is None: + iscale.set_scaling_factor(self.voltage, 1) + if hasattr(self, "conc_water") and ( + iscale.get_scaling_factor(self.conc_water, warning=True) is None + ): + iscale.set_scaling_factor(self.conc_water, 1e-4) + + if hasattr(self, "acid_produced") and ( + iscale.get_scaling_factor(self.acid_produced, warning=True) is None + ): + sf = 1 + for j in self.config.property_package.anion_set: + sf = smooth_min( + iscale.get_scaling_factor( + self.acidic.properties_in[0].flow_mass_phase_comp["Liq", j] + ), + iscale.get_scaling_factor( + self.acidic.properties_in[0].flow_mass_phase_comp["Liq", "H_+"] + ), + ) + + iscale.set_scaling_factor(self.acid_produced, sf) + + if hasattr(self, "base_produced") and ( + iscale.get_scaling_factor(self.base_produced, warning=True) is None + ): + sf = 1 + for j in self.config.property_package.cation_set: + sf = smooth_min( + iscale.get_scaling_factor( + self.basic.properties_in[0].flow_mass_phase_comp["Liq", j] + ), + iscale.get_scaling_factor( + self.basic.properties_in[0].flow_mass_phase_comp["Liq", "OH_-"] + ), + ) + + iscale.set_scaling_factor(self.base_produced, sf) + + if self.config.has_catalyst == True: + if ( + iscale.get_scaling_factor( + self.membrane_fixed_catalyst_cem, warning=True + ) + is None + ): + iscale.set_scaling_factor(self.membrane_fixed_catalyst_cem, 1e-3) + # if self.config.has_catalyst == True: + if ( + iscale.get_scaling_factor( + self.membrane_fixed_catalyst_aem, warning=True + ) + is None + ): + iscale.set_scaling_factor(self.membrane_fixed_catalyst_aem, 1e-3) + + if iscale.get_scaling_factor(self.k_a, warning=True) is None: + iscale.set_scaling_factor(self.k_a, 1e6) + + if iscale.get_scaling_factor(self.k_b, warning=True) is None: + iscale.set_scaling_factor(self.k_b, 1e2) + + if ( + self.config.has_catalyst == True + or self.config.limiting_potential_method_bpem + == LimitingpotentialMethod.Empirical + ) and iscale.get_scaling_factor(self.elec_field_non_dim, warning=True) is None: + iscale.set_scaling_factor(self.elec_field_non_dim, 1e-1) + + if ( + self.config.limiting_potential_method_bpem + == LimitingpotentialMethod.Empirical + and iscale.get_scaling_factor(self.potential_barrier_bpem, warning=True) + is None + ): + iscale.set_scaling_factor(self.potential_barrier_bpem, 1) + + if hasattr(self, "membrane_fixed_charge") and ( + iscale.get_scaling_factor(self.membrane_fixed_charge, warning=True) is None + ): + iscale.set_scaling_factor(self.membrane_fixed_charge, 1e-3) + if hasattr(self, "diffus_mass") and ( + iscale.get_scaling_factor(self.diffus_mass, warning=True) is None + ): + iscale.set_scaling_factor(self.diffus_mass, 1e9) + if hasattr(self, "salt_conc_aem") and ( + iscale.get_scaling_factor(self.salt_conc_aem, warning=True) is None + ): + iscale.set_scaling_factor(self.salt_conc_aem, 1e-2) + if hasattr(self, "salt_conc_cem") and ( + iscale.get_scaling_factor(self.salt_conc_cem, warning=True) is None + ): + iscale.set_scaling_factor(self.salt_conc_cem, 1e-2) + + if ( + self.config.has_catalyst == True + or self.config.limiting_potential_method_bpem + == LimitingpotentialMethod.Empirical + ): + + if hasattr(self, "relative_permittivity") and ( + iscale.get_scaling_factor(self.relative_permittivity, warning=True) + is None + ): + iscale.set_scaling_factor(self.relative_permittivity, 1e-1) + + if iscale.get_scaling_factor(self.kr, warning=True) is None: + iscale.set_scaling_factor(self.kr, 1e-11) + if ( + hasattr(self, "k2_zero") + and iscale.get_scaling_factor(self.k2_zero, warning=True) is None + ): + iscale.set_scaling_factor(self.k2_zero, 1e5) + + # The folloing Vars are built for constructing constraints and their sf are computed from other Vars. + + if ( + self.config.has_catalyst == True + and iscale.get_scaling_factor(self.potential_membrane_bpem, warning=True) + is None + ): + sf = ( + ( + iscale.get_scaling_factor(self.elec_field_non_dim) + * iscale.get_scaling_factor(self.relative_permittivity) + * 293**-2 + / 0.09636**-1 + ) + ** 2 + * value(Constants.vacuum_electric_permittivity) ** -1 + * iscale.get_scaling_factor(self.relative_permittivity) + * value(Constants.faraday_constant) ** -1 + * iscale.get_scaling_factor(self.membrane_fixed_charge) + ) + + iscale.set_scaling_factor(self.potential_membrane_bpem, 1e1) + + if self.config.has_catalyst == True: + if iscale.get_scaling_factor(self.flux_splitting, warning=True) is None: + + terms = 40 + sf = 0 + for indx in range(terms): + sf += ( + 2**indx + * iscale.get_scaling_factor(self.elec_field_non_dim) ** -indx + / (math.factorial(indx) * math.factorial(indx + 1)) + ) + + sf **= -1 + sf *= iscale.get_scaling_factor( + self.k2_zero + ) * iscale.get_scaling_factor(self.conc_water) + sf_a = ( + sf + * iscale.get_scaling_factor(self.membrane_fixed_catalyst_cem) + / iscale.get_scaling_factor(self.k_a) + ) + sf_b = ( + sf + * iscale.get_scaling_factor(self.membrane_fixed_catalyst_aem) + / iscale.get_scaling_factor(self.k_b) + ) + + sf = (sf_a**-1 + sf_b**-1) ** -1 * sqrt( + iscale.get_scaling_factor(self.potential_membrane_bpem) + * value(Constants.vacuum_electric_permittivity) ** -1 + * iscale.get_scaling_factor(self.relative_permittivity) + / ( + value(Constants.faraday_constant) ** -1 + * iscale.get_scaling_factor(self.membrane_fixed_charge) + ) + ) + + iscale.set_scaling_factor(self.flux_splitting, sf) + + iscale.set_scaling_factor( + self.elec_migration_bpem_flux_in, + iscale.get_scaling_factor(self.current) + * iscale.get_scaling_factor(self.cell_length) ** -1 + * iscale.get_scaling_factor(self.cell_width) ** -1 + * iscale.get_scaling_factor(self.shadow_factor) ** -1 + * 1e5, + ) + iscale.set_scaling_factor( + self.elec_migration_bpem_flux_out, + iscale.get_scaling_factor(self.current) + * iscale.get_scaling_factor(self.cell_length) ** -1 + * iscale.get_scaling_factor(self.cell_width) ** -1 + * iscale.get_scaling_factor(self.shadow_factor) ** -1 + * 1e5, + ) + + for ind in self.generation_cem_flux_in: + if ind[2] == "H_+" or "H2O": + if self.config.has_catalyst == True: + sf = 0.5 * iscale.get_scaling_factor(self.flux_splitting) + else: + sf = iscale.get_scaling_factor(self.elec_migration_bpem_flux_in) + else: + sf = 1 + + iscale.set_scaling_factor(self.generation_cem_flux_in[ind], sf) + + for ind in self.generation_cem_flux_out: + if ind[2] == "H_+" or "H2O": + if self.config.has_catalyst == True: + sf = 0.5 * iscale.get_scaling_factor(self.flux_splitting) + else: + sf = iscale.get_scaling_factor(self.elec_migration_bpem_flux_out) + else: + sf = 1 + + iscale.set_scaling_factor(self.generation_cem_flux_out[ind], sf) + + for ind in self.generation_aem_flux_in: + if ind[2] == "OH_-" or "H2O": + if self.config.has_catalyst == True: + sf = iscale.get_scaling_factor(self.flux_splitting) + else: + sf = iscale.get_scaling_factor(self.elec_migration_bpem_flux_in) + else: + sf = 1 + + iscale.set_scaling_factor(self.generation_aem_flux_in[ind], sf) + + for ind in self.generation_aem_flux_out: + if ind[2] == "OH_-" or "H2O": + if self.config.has_catalyst == True: + sf = iscale.get_scaling_factor(self.flux_splitting) + else: + sf = iscale.get_scaling_factor(self.elec_migration_bpem_flux_out) + else: + sf = 1 + + iscale.set_scaling_factor(self.generation_aem_flux_out[ind], sf) + + for ind in self.nonelec_bpem_flux_in: + if ind[2] == "H2O": + sf = ( + 1e-3 + * 0.018 + * iscale.get_scaling_factor(self.water_permeability_membrane) + * iscale.get_scaling_factor( + self.acidic.properties_in[ind[0]].pressure_osm_phase[ind[1]] + ) + ) + else: + sf = 1 + iscale.set_scaling_factor(self.nonelec_bpem_flux_in[ind], sf) + for ind in self.nonelec_bpem_flux_out: + if ind[2] == "H2O": + sf = ( + 1e-3 + * 0.018 + * iscale.get_scaling_factor(self.water_permeability_membrane) + * iscale.get_scaling_factor( + self.acidic.properties_out[ind[0]].pressure_osm_phase[ind[1]] + ) + ) + else: + sf = 1 + + iscale.set_scaling_factor(self.nonelec_bpem_flux_out[ind], sf) + + for ind in self.acidic.mass_transfer_term: + if ind[2] == "H_+": + sf = iscale.get_scaling_factor(self.generation_cem_flux_in[ind]) + else: + if ind[2] == "H2O": + sf = iscale.get_scaling_factor(self.nonelec_bpem_flux_in[ind]) + else: + sf = iscale.get_scaling_factor( + self.elec_migration_bpem_flux_in[ind] + ) + + sf *= ( + iscale.get_scaling_factor(self.cell_width) + * iscale.get_scaling_factor(self.shadow_factor) + * iscale.get_scaling_factor(self.cell_length) + * iscale.get_scaling_factor(self.cell_num) + ) + iscale.set_scaling_factor(self.acidic.mass_transfer_term[ind], sf) + # + for ind in self.basic.mass_transfer_term: + if ind[2] == "OH_-": + sf = iscale.get_scaling_factor(self.generation_aem_flux_in[ind]) + else: + if ind[2] == "H2O": + sf = iscale.get_scaling_factor(self.nonelec_bpem_flux_in[ind]) + else: + sf = iscale.get_scaling_factor( + self.elec_migration_bpem_flux_in[ind] + ) + + sf *= ( + iscale.get_scaling_factor(self.cell_width) + * iscale.get_scaling_factor(self.shadow_factor) + * iscale.get_scaling_factor(self.cell_length) + * iscale.get_scaling_factor(self.cell_num) + ) + iscale.set_scaling_factor(self.basic.mass_transfer_term[ind], sf) + + for ind, c in self.specific_power_electrical.items(): + iscale.set_scaling_factor( + self.specific_power_electrical[ind], + 3.6e6 + * iscale.get_scaling_factor(self.current[ind]) + * iscale.get_scaling_factor(self.voltage[ind]) + * iscale.get_scaling_factor( + self.acidic.properties_out[ind].flow_vol_phase["Liq"] + ) + ** -1, + ) + + for ind in self.velocity_basic: + if ( + iscale.get_scaling_factor(self.velocity_basic[ind], warning=False) + is None + ): + sf = ( + iscale.get_scaling_factor( + self.basic.properties_in[ind].flow_vol_phase["Liq"] + ) + * iscale.get_scaling_factor(self.cell_width) ** -1 + * iscale.get_scaling_factor(self.shadow_factor) ** -1 + * iscale.get_scaling_factor(self.channel_height) ** -1 + * iscale.get_scaling_factor(self.spacer_porosity) ** -1 + * iscale.get_scaling_factor(self.cell_num) ** -1 + ) + + iscale.set_scaling_factor(self.velocity_basic[ind], sf) + + for ind in self.velocity_acidic: + if ( + iscale.get_scaling_factor(self.velocity_acidic[ind], warning=False) + is None + ): + sf = ( + iscale.get_scaling_factor( + self.acidic.properties_in[ind].flow_vol_phase["Liq"] + ) + * iscale.get_scaling_factor(self.cell_width) ** -1 + * iscale.get_scaling_factor(self.shadow_factor) ** -1 + * iscale.get_scaling_factor(self.channel_height) ** -1 + * iscale.get_scaling_factor(self.spacer_porosity) ** -1 + * iscale.get_scaling_factor(self.cell_num) ** -1 + ) + + iscale.set_scaling_factor(self.velocity_acidic[ind], sf) + + if hasattr(self, "spacer_specific_area") and ( + iscale.get_scaling_factor(self.spacer_specific_area, warning=True) is None + ): + iscale.set_scaling_factor(self.spacer_specific_area, 1e-4) + if hasattr(self, "hydraulic_diameter") and ( + iscale.get_scaling_factor(self.hydraulic_diameter, warning=True) is None + ): + iscale.set_scaling_factor(self.hydraulic_diameter, 1e4) + if hasattr(self, "N_Re") and ( + iscale.get_scaling_factor(self.N_Re, warning=True) is None + ): + sf = ( + iscale.get_scaling_factor(self.dens_mass) + * iscale.get_scaling_factor(self.velocity_acidic[0]) + * iscale.get_scaling_factor(self.hydraulic_diameter) + * iscale.get_scaling_factor(self.visc_d) ** -1 + ) + iscale.set_scaling_factor(self.N_Re, sf) + if hasattr(self, "N_Sc") and ( + iscale.get_scaling_factor(self.N_Sc, warning=True) is None + ): + sf = ( + iscale.get_scaling_factor(self.visc_d) + * iscale.get_scaling_factor(self.dens_mass) ** -1 + * iscale.get_scaling_factor(self.diffus_mass) ** -1 + ) + iscale.set_scaling_factor(self.N_Sc, sf) + if hasattr(self, "N_Sh") and ( + iscale.get_scaling_factor(self.N_Sh, warning=True) is None + ): + sf = ( + 10 + * iscale.get_scaling_factor(self.N_Re) ** 0.5 + * iscale.get_scaling_factor(self.N_Sc) ** 0.33 + ) + iscale.set_scaling_factor(self.N_Sh, sf) + if hasattr(self, "friction_factor") and ( + iscale.get_scaling_factor(self.friction_factor, warning=True) is None + ): + if self.config.friction_factor_method == FrictionFactorMethod.fixed: + sf = 0.1 + elif self.config.friction_factor_method == FrictionFactorMethod.Gurreri: + sf = ( + (4 * 50.6) ** -1 + * (iscale.get_scaling_factor(self.spacer_porosity)) ** -7.06 + * iscale.get_scaling_factor(self.N_Re) ** -1 + ) + elif self.config.friction_factor_method == FrictionFactorMethod.Kuroda: + sf = (4 * 9.6) ** -1 * iscale.get_scaling_factor(self.N_Re) ** -0.5 + iscale.set_scaling_factor(self.friction_factor, sf) + + if hasattr(self, "pressure_drop") and ( + iscale.get_scaling_factor(self.pressure_drop, warning=True) is None + ): + if self.config.pressure_drop_method == PressureDropMethod.experimental: + sf = 1e-5 + else: + sf = ( + iscale.get_scaling_factor(self.dens_mass) + * iscale.get_scaling_factor(self.friction_factor) + * iscale.get_scaling_factor(self.velocity_acidic[0]) ** 2 + * 2 + * iscale.get_scaling_factor(self.hydraulic_diameter) ** -1 + ) + iscale.set_scaling_factor(self.pressure_drop, sf) + + if hasattr(self, "pressure_drop_total") and ( + iscale.get_scaling_factor(self.pressure_drop_total, warning=True) is None + ): + if self.config.pressure_drop_method == PressureDropMethod.experimental: + sf = 1e-5 * iscale.get_scaling_factor(self.cell_length) + else: + sf = ( + iscale.get_scaling_factor(self.dens_mass) + * iscale.get_scaling_factor(self.friction_factor) + * iscale.get_scaling_factor(self.velocity_acidic[0]) ** 2 + * 2 + * iscale.get_scaling_factor(self.hydraulic_diameter) ** -1 + * iscale.get_scaling_factor(self.cell_length) + ) + iscale.set_scaling_factor(self.pressure_drop_total, sf) + + if hasattr(self, "current_dens_lim_bpem"): + if iscale.get_scaling_factor(self.current_dens_lim_bpem) is None: + if ( + self.config.limiting_current_density_method_bpem + == LimitingCurrentDensitybpemMethod.InitialValue + ): + sf = self.config.limiting_current_density_bpem_data**-1 + iscale.set_scaling_factor(self.current_dens_lim_bpem, sf) + elif ( + self.config.limiting_current_density_method_bpem + == LimitingCurrentDensitybpemMethod.Empirical + ): + sf = ( + iscale.get_scaling_factor(self.diffus_mass) + * value(Constants.faraday_constant) ** -1 + * (1 / 2 * iscale.get_scaling_factor(self.salt_conc_aem)) ** 2 + / ( + iscale.get_scaling_factor(self.membrane_thickness) + * iscale.get_scaling_factor(self.membrane_fixed_charge) + ) + ) + + iscale.set_scaling_factor(self.current_dens_lim_bpem, sf) + + iscale.set_scaling_factor( + self.power_electrical, + iscale.get_scaling_factor(self.current) + * iscale.get_scaling_factor(self.voltage), + ) + + # Constraint scaling + for ind, c in self.eq_current_voltage_relation.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.membrane_areal_resistance) + ) + for ind, c in self.eq_power_electrical.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.power_electrical) + ) + for ind, c in self.eq_specific_power_electrical.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.specific_power_electrical[ind]) + ) + for ind, c in self.eq_elec_migration_bpem_flux_in.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.elec_migration_bpem_flux_in) + ) + for ind, c in self.eq_elec_migration_bpem_flux_out.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.elec_migration_bpem_flux_in) + ) + + for ind, c in self.eq_nonelec_bpem_flux_in.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.nonelec_bpem_flux_in[ind]) + ) + + for ind, c in self.eq_nonelec_bpem_flux_out.items(): + iscale.constraint_scaling_transform( + c, iscale.get_scaling_factor(self.nonelec_bpem_flux_out[ind]) + ) + + for ind, c in self.eq_mass_transfer_term_basic.items(): + iscale.constraint_scaling_transform( + c, + min( + iscale.get_scaling_factor(self.generation_aem_flux_in[ind]), + iscale.get_scaling_factor(self.elec_migration_bpem_flux_in[ind]), + iscale.get_scaling_factor( + self.nonelec_bpem_flux_in[ind], + self.elec_migration_bpem_flux_out[ind], + ), + iscale.get_scaling_factor(self.nonelec_bpem_flux_out[ind]), + ) + * iscale.get_scaling_factor(self.cell_width) + * iscale.get_scaling_factor(self.shadow_factor) + * iscale.get_scaling_factor(self.cell_length) + * iscale.get_scaling_factor(self.cell_num), + ) + for ind, c in self.eq_mass_transfer_term_acidic.items(): + iscale.constraint_scaling_transform( + c, + min( + iscale.get_scaling_factor(self.generation_cem_flux_in[ind]), + iscale.get_scaling_factor(self.elec_migration_bpem_flux_out[ind]), + iscale.get_scaling_factor( + self.nonelec_bpem_flux_in[ind], + self.elec_migration_bpem_flux_out[ind], + ), + iscale.get_scaling_factor(self.nonelec_bpem_flux_out[ind]), + ) + * iscale.get_scaling_factor(self.cell_width) + * iscale.get_scaling_factor(self.shadow_factor) + * iscale.get_scaling_factor(self.cell_length) + * iscale.get_scaling_factor(self.cell_num), + ) + + for ind, c in self.eq_power_electrical.items(): + iscale.constraint_scaling_transform( + c, + iscale.get_scaling_factor(self.power_electrical[ind]), + ) + + for ind, c in self.eq_specific_power_electrical.items(): + iscale.constraint_scaling_transform( + c, + iscale.get_scaling_factor(self.specific_power_electrical[ind]) + * iscale.get_scaling_factor( + self.acidic.properties_out[ind].flow_vol_phase["Liq"] + ), + ) + + def _get_stream_table_contents(self, time_point=0): + return create_stream_table_dataframe( + { + "base channel of the bipolar membrane Channel Inlet": self.inlet_basic, + "acid channel of the bipolar membrane Channel Inlet": self.inlet_acidic, + "base channel of the bipolar membrane Channel Outlet": self.outlet_basic, + "acid channel of the bipolar membrane Channel Outlet": self.outlet_acidic, + }, + time_point=time_point, + ) + + def _get_performance_contents(self, time_point=0): + return { + "vars": { + "Total electrical power consumption(Watt)": self.power_electrical[ + time_point + ], + "Specific electrical power consumption (kW*h/m**3)": self.specific_power_electrical[ + time_point + ], + }, + "exprs": {}, + "params": {}, + } + + def get_power_electrical(self, time_point=0): + return self.power_electrical[time_point] + + @property + def default_costing_method(self): + return cost_bipolar_electrodialysis diff --git a/watertap/unit_models/tests/test_biploar_electrodialysis_0D.py b/watertap/unit_models/tests/test_biploar_electrodialysis_0D.py new file mode 100644 index 0000000000..37d90c2f45 --- /dev/null +++ b/watertap/unit_models/tests/test_biploar_electrodialysis_0D.py @@ -0,0 +1,746 @@ +################################################################################# +# WaterTAP Copyright (c) 2020-2024, The Regents of the University of California, +# through Lawrence Berkeley National Laboratory, Oak Ridge National Laboratory, +# National Renewable Energy Laboratory, and National Energy Technology +# Laboratory (subject to receipt of any required approvals from the U.S. Dept. +# of Energy). All rights reserved. +# +# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license +# information, respectively. These files are also available online at the URL +# "https://github.com/watertap-org/watertap/" +################################################################################# +import numpy as np +import idaes.core.util.scaling as iscale +import pytest +from idaes.core import ( + FlowsheetBlock, +) + +from idaes.core.util.model_statistics import degrees_of_freedom +from idaes.core.util.testing import initialization_tester +from pyomo.environ import ( + ConcreteModel, + assert_optimal_termination, + value, +) + +from watertap.unit_models.Bipolar_Electrodialysis_0D import ( + Bipolar_Electrodialysis_0D, + LimitingCurrentDensitybpemMethod, + LimitingpotentialMethod, + PressureDropMethod, + FrictionFactorMethod, + HydraulicDiameterMethod, +) +from watertap.core.solvers import get_solver +from watertap.property_models.multicomp_aq_sol_prop_pack import MCASParameterBlock + + +__author__ = "Johnson Dhanasekaran" + +solver = get_solver() + + +# ----------------------------------------------------------------------------- +# Start test class + + +class Test_catalyst: + @pytest.fixture(scope="class") + def bped(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + has_catalyst=True, + limiting_current_density_method_bpem=LimitingCurrentDensitybpemMethod.Empirical, + ) + + m.fs.unit.diffus_mass.fix((2.03 + 1.96) * 10**-9 / 2) + m.fs.unit.membrane_fixed_charge["bpem"].fix(5e3) + m.fs.unit.salt_conc_aem["bpem"].fix(2000) + m.fs.unit.salt_conc_cem["bpem"].fix(2000) + m.fs.unit.conc_water["bpem"].fix(50 * 1e3) + m.fs.unit.kr["bpem"].fix(1.3 * 10**10) + m.fs.unit.k2_zero["bpem"].fix(2 * 10**-6) + m.fs.unit.relative_permittivity["bpem"].fix(30) + m.fs.unit.membrane_fixed_catalyst_cem["bpem"].fix(5e3) + m.fs.unit.membrane_fixed_catalyst_aem["bpem"].fix(5e3) + m.fs.unit.k_a.fix(447) + m.fs.unit.k_b.fix(5e4) + + m.fs.unit.ion_trans_number_membrane["bpem", "Na_+"].fix(0.5) + m.fs.unit.ion_trans_number_membrane["bpem", "Cl_-"].fix(0.5) + m.fs.unit.ion_trans_number_membrane["bpem", "H_+"].fix(0.1) + m.fs.unit.ion_trans_number_membrane["bpem", "OH_-"].fix(0.1) + + # Set inlet streams. + m.fs.unit.inlet_basic.pressure.fix(101325) + m.fs.unit.inlet_basic.temperature.fix(298.15) + m.fs.unit.inlet_acidic.pressure.fix(101325) + m.fs.unit.inlet_acidic.temperature.fix(298.15) + m.fs.unit.spacer_porosity.fix(1) + + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-4) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-4) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "H_+"].fix(7.38e-4) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "OH_-"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "H_+"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "OH_-"].fix(7.38e-4) + + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + + m.fs.unit.shadow_factor.fix(1) + m.fs.unit.water_trans_number_membrane["bpem"].fix((5.8 + 4.3) / 2) + m.fs.unit.water_permeability_membrane["bpem"].fix((2.16e-14 + 1.75e-14) / 2) + m.fs.unit.electrodes_resistance.fix(0) + m.fs.unit.current_utilization.fix(1) + m.fs.unit.channel_height.fix(2.7e-4) + m.fs.unit.membrane_areal_resistance.fix((1.89e-4 + 1.77e-4) / 2) + m.fs.unit.cell_width.fix(0.1) + m.fs.unit.cell_length.fix(0.79) + m.fs.unit.membrane_thickness["bpem"].fix(8e-4) + + # Set scaling of critical quantities. + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e3, index=("Liq", "Na_+") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e3, index=("Liq", "Cl_-") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e3, index=("Liq", "H_+") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e3, index=("Liq", "OH_-") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e2, index=("Liq", "H2O") + ) + m.fs.unit.cell_num.fix(1) + iscale.set_scaling_factor(m.fs.unit.k_a, 1e-2) + iscale.set_scaling_factor(m.fs.unit.k_b, 1e-4) + iscale.set_scaling_factor(m.fs.unit.voltage, 1e-1) + iscale.set_scaling_factor(m.fs.unit.flux_splitting, 1e3) + return m + + @pytest.mark.unit + def test_assign(self, bped): + + # Experimental data from Wilhelm et al. (2002) with additional inputs from Mareev et al. (2020) + expt_membrane_potential = np.array( + [0.088, 0.184, 0.4045, 0.690, 0.858, 0.95, 1.3] + ) # in volts + + # experimental current density: 11.7, 15.5, 20.8, 24.4, 30.6, 40.0, 100.6] # in mA/cm2 + expected_current_density = np.array( + [19.250, 19.269, 19.606, 22.764, 29.009, 35.278, 99.223] + ) # in mA/cm2 (computed numerically) + + for indx, v in enumerate(expt_membrane_potential): + m = bped + m.fs.unit.potential_membrane_bpem[0].fix(v) + + iscale.calculate_scaling_factors(m.fs) + assert degrees_of_freedom(m) == 0 + initialization_tester(m) + results = solver.solve(m) + assert_optimal_termination(results) + + current_density_computed = ( + 0.1 + * m.fs.unit.current[0] + / (m.fs.unit.cell_width * m.fs.unit.cell_length) + ) # convert to mA/cm2 + assert value(current_density_computed) == pytest.approx( + expected_current_density[indx], rel=1e-3 + ) + + +class Test_limiting_parameters: + @pytest.fixture(scope="class") + def limiting_current_check(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + limiting_current_density_method_bpem=LimitingCurrentDensitybpemMethod.Empirical, + limiting_potential_method_bpem=LimitingpotentialMethod.InitialValue, + limiting_potential_data=0.5, + has_catalyst=False, + ) + return m + + @pytest.fixture(scope="class") + def potential_barrier_check(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + limiting_current_density_method_bpem=LimitingCurrentDensitybpemMethod.Empirical, + limiting_potential_method_bpem=LimitingpotentialMethod.Empirical, + has_catalyst=False, + ) + + m.fs.unit.kr["bpem"].fix(1.33 * 10**11) + m.fs.unit.k2_zero["bpem"].fix(2 * 10**-5) + m.fs.unit.relative_permittivity["bpem"].fix(20) + + return m + + @pytest.mark.unit + def test_assign(self, limiting_current_check, potential_barrier_check): + check_m = (limiting_current_check, potential_barrier_check) + for m in check_m: + m.fs.unit.shadow_factor.fix(1) + m.fs.unit.current.fix(1e1) + m.fs.unit.water_trans_number_membrane["bpem"].fix((5.8 + 4.3) / 2) + m.fs.unit.water_permeability_membrane["bpem"].fix((2.16e-14 + 1.75e-14) / 2) + m.fs.unit.electrodes_resistance.fix(0) + m.fs.unit.cell_num.fix(1) + m.fs.unit.current_utilization.fix(1) + m.fs.unit.channel_height.fix(2.7e-4) + m.fs.unit.membrane_areal_resistance.fix((1.89e-4 + 1.77e-4) / 2) + m.fs.unit.cell_width.fix(0.1) + m.fs.unit.cell_length.fix(0.79) + m.fs.unit.membrane_thickness["bpem"].fix(1.3e-4) + m.fs.unit.diffus_mass.fix((2.03 + 1.96) * 10**-9 / 2) + m.fs.unit.ion_trans_number_membrane["bpem", "Na_+"].fix(0.5) + m.fs.unit.ion_trans_number_membrane["bpem", "Cl_-"].fix(0.5) + m.fs.unit.ion_trans_number_membrane["bpem", "H_+"].fix(0.1) + m.fs.unit.ion_trans_number_membrane["bpem", "OH_-"].fix(0.1) + m.fs.unit.conc_water["bpem"].fix(55 * 1e3) + m.fs.unit.kr["bpem"].fix(1.33 * 10**11) + m.fs.unit.k2_zero["bpem"].fix(2 * 10**-5) + m.fs.unit.relative_permittivity["bpem"].fix(20) + m.fs.unit.diffus_mass.fix((2.03 + 1.96) * 10**-9 / 2) + m.fs.unit.salt_conc_aem["bpem"].fix(500 + 2 * 250) + m.fs.unit.salt_conc_cem["bpem"].fix(500 + 2 * 250) + m.fs.unit.membrane_fixed_charge["bpem"].fix(1.5e3) + + # Set inlet streams. + m.fs.unit.inlet_basic.pressure.fix(101325) + m.fs.unit.inlet_basic.temperature.fix(298.15) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-1) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-1) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "H_+"].fix(7.38e-1) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "OH_-"].fix(7.38e-1) + m.fs.unit.inlet_acidic.pressure.fix(101325) + m.fs.unit.inlet_acidic.temperature.fix(298.15) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-1) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-1) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "H_+"].fix(7.38e-1) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "OH_-"].fix(7.38e-1) + m.fs.unit.spacer_porosity.fix(1) + + # Set scaling of critical quantities. + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e1, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "Na_+") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "Cl_-") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "H_+") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "OH_-") + ) + + # Data on limiting current and potential barrier to water splitting have been obtained from: + # Fumatech, Technical Data Sheet for Fumasep FBM, 2020. With additional modelling parameters obtained from + # Ionescu, Viorel. Advanced Topics in Optoelectronics, Microelectronics, and Nanotechnologies (2023) + # Expected current density is 1000 A/m2 and potential barrier is 0.8 V + + iscale.calculate_scaling_factors(check_m[1]) + + # Test computing limiting current in bipolar membrane + iscale.calculate_scaling_factors(check_m[0]) + assert degrees_of_freedom(check_m[0]) == 0 + initialization_tester(check_m[0]) + results = solver.solve(check_m[0]) + assert_optimal_termination(results) + assert value(check_m[0].fs.unit.current_dens_lim_bpem[0]) == pytest.approx( + 987.120, rel=1e-3 + ) + + # Test computing limiting current and potential barrier to water splitting in bipolar membrane + iscale.calculate_scaling_factors(check_m[1]) + assert degrees_of_freedom(check_m[1]) == 0 + initialization_tester(check_m[1]) + results = solver.solve(check_m[1]) + assert_optimal_termination(results) + assert value(check_m[0].fs.unit.current_dens_lim_bpem[0]) == pytest.approx( + 987.120, rel=1e-3 + ) + assert value(check_m[1].fs.unit.potential_barrier_bpem[0]) == pytest.approx( + 0.786, rel=1e-3 + ) + + +class Test_BPED_pressure_drop_components: + + @pytest.fixture(scope="class") + def bped_m0(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + has_catalyst=False, + pressure_drop_method=PressureDropMethod.experimental, + has_pressure_change=True, + ) + return m + + @pytest.fixture(scope="class") + def bped_m1(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + has_catalyst=False, + pressure_drop_method=PressureDropMethod.Darcy_Weisbach, + friction_factor_method=FrictionFactorMethod.fixed, + hydraulic_diameter_method=HydraulicDiameterMethod.conventional, + has_pressure_change=True, + ) + return m + + @pytest.fixture(scope="class") + def bped_m2(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + pressure_drop_method=PressureDropMethod.Darcy_Weisbach, + friction_factor_method=FrictionFactorMethod.Gurreri, + hydraulic_diameter_method=HydraulicDiameterMethod.conventional, + has_pressure_change=True, + ) + return m + + @pytest.fixture(scope="class") + def bped_m3(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + pressure_drop_method=PressureDropMethod.Darcy_Weisbach, + friction_factor_method=FrictionFactorMethod.Kuroda, + hydraulic_diameter_method=HydraulicDiameterMethod.conventional, + has_pressure_change=True, + ) + return m + + @pytest.fixture(scope="class") + def bped_m4(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + pressure_drop_method=PressureDropMethod.Darcy_Weisbach, + friction_factor_method=FrictionFactorMethod.Kuroda, + hydraulic_diameter_method=HydraulicDiameterMethod.fixed, + has_pressure_change=True, + ) + return m + + @pytest.fixture(scope="class") + def bped_m5(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-", "H_+", "OH_-"], + "mw_data": { + "Na_+": 23e-3, + "Cl_-": 35.5e-3, + "H_+": 1e-3, + "OH_-": 17.0e-3, + }, + "elec_mobility_data": { + ("Liq", "Na_+"): 5.19e-8, + ("Liq", "Cl_-"): 7.92e-8, + ("Liq", "H_+"): 36.23e-8, + ("Liq", "OH_-"): 20.64e-8, + }, + "charge": {"Na_+": 1, "Cl_-": -1, "H_+": 1, "OH_-": -1}, + "diffusivity_data": { + ("Liq", "Na_+"): 1.33e-9, + ("Liq", "Cl_-"): 2.03e-9, + ("Liq", "H_+"): 9.31e-9, + ("Liq", "OH_-"): 5.27e-9, + }, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Bipolar_Electrodialysis_0D( + property_package=m.fs.properties, + pressure_drop_method=PressureDropMethod.Darcy_Weisbach, + friction_factor_method=FrictionFactorMethod.Kuroda, + hydraulic_diameter_method=HydraulicDiameterMethod.spacer_specific_area_known, + has_pressure_change=True, + ) + return m + + @pytest.mark.unit + def test_deltaP_various_methods( + self, bped_m0, bped_m1, bped_m2, bped_m3, bped_m4, bped_m5 + ): + bped_m = (bped_m0, bped_m1, bped_m2, bped_m3, bped_m4, bped_m5) + for m in bped_m: + m.fs.unit.shadow_factor.fix(1) + m.fs.unit.current.fix(1e2) + m.fs.unit.water_trans_number_membrane["bpem"].fix((5.8 + 4.3) / 2) + m.fs.unit.water_permeability_membrane["bpem"].fix((2.16e-14 + 1.75e-14) / 2) + m.fs.unit.electrodes_resistance.fix(0) + m.fs.unit.cell_num.fix(10) + m.fs.unit.current_utilization.fix(1) + m.fs.unit.channel_height.fix(2.7e-4) + m.fs.unit.membrane_areal_resistance.fix((1.89e-4 + 1.77e-4) / 2) + m.fs.unit.cell_width.fix(0.1) + m.fs.unit.cell_length.fix(0.79) + m.fs.unit.membrane_thickness["bpem"].fix(1.3e-4) + m.fs.unit.diffus_mass.fix((2.03 + 1.96) * 10**-9 / 2) + m.fs.unit.ion_trans_number_membrane["bpem", "Na_+"].fix(0.5) + m.fs.unit.ion_trans_number_membrane["bpem", "Cl_-"].fix(0.5) + m.fs.unit.ion_trans_number_membrane["bpem", "H_+"].fix(0.1) + m.fs.unit.ion_trans_number_membrane["bpem", "OH_-"].fix(0.1) + m.fs.unit.conc_water["bpem"].fix(55 * 1e3) + m.fs.unit.kr["bpem"].fix(1.33 * 10**11) + m.fs.unit.k2_zero["bpem"].fix(2 * 10**-5) + m.fs.unit.relative_permittivity["bpem"].fix(20) + m.fs.unit.diffus_mass.fix((2.03 + 1.96) * 10**-9 / 2) + m.fs.unit.salt_conc_aem["bpem"].fix(500 + 2 * 250) + m.fs.unit.salt_conc_cem["bpem"].fix(500 + 2 * 250) + m.fs.unit.membrane_fixed_charge["bpem"].fix(1.5e3) + + # Set inlet streams. + m.fs.unit.inlet_basic.pressure.fix(201035) + m.fs.unit.inlet_basic.temperature.fix(298.15) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-4) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-4) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "H_+"].fix(7.38e-4) + m.fs.unit.inlet_basic.flow_mol_phase_comp[0, "Liq", "OH_-"].fix(7.38e-4) + m.fs.unit.inlet_acidic.pressure.fix(201035) + m.fs.unit.inlet_acidic.temperature.fix(298.15) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "H_+"].fix(7.38e-4) + m.fs.unit.inlet_acidic.flow_mol_phase_comp[0, "Liq", "OH_-"].fix(7.38e-4) + m.fs.unit.spacer_porosity.fix(0.83) + + # Set scaling of critical quantities. + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "H2O") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "Na_+") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "Cl_-") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "H_+") + ) + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e0, index=("Liq", "OH_-") + ) + + iscale.set_scaling_factor(m.fs.unit.cell_length, 1) + iscale.set_scaling_factor(m.fs.unit.cell_num, 0.1) + iscale.set_scaling_factor(m.fs.unit.voltage, 1e-1) + + # Test bped_m0 + bped_m[0].fs.unit.pressure_drop.fix(1e4) + iscale.calculate_scaling_factors(bped_m[0]) + assert degrees_of_freedom(bped_m[0]) == 0 + initialization_tester(bped_m[0]) + results = solver.solve(bped_m[0]) + assert_optimal_termination(results) + assert value(bped_m[0].fs.unit.pressure_drop_total[0]) == pytest.approx( + 7900, rel=1e-3 + ) + + # Test bped_m1 + bped_m[1].fs.unit.current.fix(1e2) + iscale.set_scaling_factor(bped_m[1].fs.unit.current, 1e-2) + bped_m[1].fs.unit.diffus_mass.fix(1.6e-9) + bped_m[1].fs.unit.friction_factor.fix(20) + iscale.calculate_scaling_factors(bped_m[1]) + assert degrees_of_freedom(bped_m[1]) == 0 + initialization_tester(bped_m[1]) + results = solver.solve(bped_m[1]) + assert_optimal_termination(results) + assert value(bped_m[1].fs.unit.N_Re) == pytest.approx(9.585, rel=1e-3) + + assert value(bped_m[1].fs.unit.pressure_drop[0]) == pytest.approx( + 10286.288, rel=1e-3 + ) + + assert value(bped_m[1].fs.unit.pressure_drop_total[0]) == pytest.approx( + 8126.168, rel=1e-3 + ) + + # Test bped_m2 + bped_m[2].fs.unit.diffus_mass.fix(1.6e-9) + iscale.calculate_scaling_factors(bped_m[2]) + assert degrees_of_freedom(bped_m[2]) == 0 + initialization_tester(bped_m[2]) + results = solver.solve(bped_m[2]) + assert_optimal_termination(results) + assert value(bped_m[2].fs.unit.N_Re) == pytest.approx(9.585, rel=1e-3) + + assert value(bped_m[2].fs.unit.pressure_drop[0]) == pytest.approx( + 40473.174, rel=1e-3 + ) + + assert value(bped_m[2].fs.unit.pressure_drop_total[0]) == pytest.approx( + 31973.808, rel=1e-3 + ) + + # Test bped_m3 + bped_m[3].fs.unit.diffus_mass.fix(1.6e-9) + iscale.calculate_scaling_factors(bped_m[3]) + assert degrees_of_freedom(bped_m[3]) == 0 + initialization_tester(bped_m[3]) + results = solver.solve(bped_m[3]) + assert_optimal_termination(results) + assert value(bped_m[3].fs.unit.N_Re) == pytest.approx(9.585, rel=1e-3) + + assert value(bped_m[3].fs.unit.pressure_drop[0]) == pytest.approx( + 7685.843, rel=1e-3 + ) + + assert value(bped_m[3].fs.unit.pressure_drop_total[0]) == pytest.approx( + 6071.816, rel=1e-3 + ) + + # Test bped_m4 + bped_m[4].fs.unit.diffus_mass.fix(1.6e-9) + bped_m[4].fs.unit.hydraulic_diameter.fix(1e-3) + iscale.calculate_scaling_factors(bped_m[4]) + assert degrees_of_freedom(bped_m[4]) == 0 + initialization_tester(bped_m[4]) + results = solver.solve(bped_m[4]) + assert_optimal_termination(results) + assert value(bped_m[4].fs.unit.N_Re) == pytest.approx(21.443, rel=1e-3) + + assert value(bped_m[4].fs.unit.pressure_drop[0]) == pytest.approx( + 2296.904, rel=1e-3 + ) + + assert value(bped_m[4].fs.unit.pressure_drop_total[0]) == pytest.approx( + 1814.554, rel=1e-3 + ) + + # Test bped_m5 + bped_m[5].fs.unit.diffus_mass.fix(1.6e-9) + bped_m[5].fs.unit.spacer_specific_area.fix(10700) + iscale.calculate_scaling_factors(bped_m[5]) + assert degrees_of_freedom(bped_m[5]) == 0 + initialization_tester(bped_m[5]) + iscale.calculate_scaling_factors(bped_m[5]) + results = solver.solve(bped_m[5]) + assert_optimal_termination(results) + assert value(bped_m[5].fs.unit.N_Re) == pytest.approx(7.716, rel=1e-3) + + assert value(bped_m[5].fs.unit.pressure_drop[0]) == pytest.approx( + 10641.053, rel=1e-3 + ) + + assert value(bped_m[5].fs.unit.pressure_drop_total[0]) == pytest.approx( + 8406.432, rel=1e-3 + )