diff --git a/docs/technical_reference/unit_models/electrodialysis_0D.rst b/docs/technical_reference/unit_models/electrodialysis_0D.rst index 5913e88ceb..fe4f2376d8 100644 --- a/docs/technical_reference/unit_models/electrodialysis_0D.rst +++ b/docs/technical_reference/unit_models/electrodialysis_0D.rst @@ -204,13 +204,16 @@ closer to the non-ideal physical conditions that can be encountered in real desa .. csv-table:: **Table 5** Essential equations supporting model extensions :header: "Description", "Equation", "Condition" - "Nonohmic potential, membrane", ":math:`\phi_m=\frac{RT}{F} \left( t_+^{iem} - t_-^{iem} \right) \ln \left( \frac{c_s^R}{c_s^L} \right)`", "`has_nonohmic_potential_membrane == True`" - "Ohmic potential, NDL", ":math:`\phi_d^{ohm}=\frac{FD_b}{\left(t_+^{iem}-t_+\right)\lambda}\ln\left(\frac{c_s^Lc_b^R}{c_s^Rc_b^L}\right)`", "`has_Nernst_diffusion_layer==True`" - "Nonohmic potential, NDL", ":math:`\phi_d^{nonohm}=\frac{RT}{F}\left(t_+-t_-\right) \ln\left(\frac{c_s^Lc_b^R}{c_s^Rc_b^L}\right)`", "`has_Nernst_diffusion_layer==True`" - "NDL thickness, cem", ":math:`\Delta^{L/R} = \frac{F D_b c_b^{L/R}}{\left(t_+^{iem}-t_+ \right) i_{lim}}`", "`has_Nernst_diffusion_layer==True`" - "NDL thickness, aem", ":math:`\Delta^{L/R} = - \frac{F D_b c_b^{L/R}}{\left(t_+^{iem}-t_+\right) i_{lim}}`", "`has_Nernst_diffusion_layer==True`" - "Concentration polarization ratio, cem", ":math:`\frac{c_s^L}{c_b^L} = 1+\frac{i}{i_{lim}},\qquad \frac{c_s^R}{c_b^R} = 1-\frac{i}{i_{lim}}`", "`has_Nernst_diffusion_layer==True` \ :sup:`1`" - "Concentration polarization ratio, aem", ":math:`\frac{c_s^L}{c_b^L} = 1-\frac{i}{i_{lim}},\qquad \frac{c_s^R}{c_b^R} = 1+\frac{i}{i_{lim}}`", "`has_Nernst_diffusion_layer==True`" + "Limiting current density", ":math:`i_{lim} = i_{lim,0}\frac{c^D}{c^D_{in}}`", "`has_Nernst_diffusion_layer==True` and `limiting_current_density_method == LimitingCurrentDensityMethod.InitialValue`" + " ", ":math:`i_{lim} = A v^B c^D`", "`has_Nernst_diffusion_layer==True` and `limiting_current_density_method == LimitingCurrentDensityMethod.Empirical`" + " ", ":math:`i_{lim} = \frac{Sh F D_b c^D}{d_H \left(t_+^{cem}-t_+ \right)}`", "`has_Nernst_diffusion_layer==True` and `limiting_current_density_method == LimitingCurrentDensityMethod.Theoretical`" + "Nonohmic potential, membrane", ":math:`\phi_m=\frac{RT}{F} \left( t_+^{iem} - t_-^{iem} \right) \ln \left( \frac{c_s^R}{c_s^L} \right)`", "`has_nonohmic_potential_membrane == True`" + "Ohmic potential, NDL", ":math:`\phi_d^{ohm}=\frac{FD_b}{\left(t_+^{iem}-t_+\right)\lambda}\ln\left(\frac{c_s^Lc_b^R}{c_s^Rc_b^L}\right)`", "`has_Nernst_diffusion_layer==True`" + "Nonohmic potential, NDL", ":math:`\phi_d^{nonohm}=\frac{RT}{F}\left(t_+-t_-\right) \ln\left(\frac{c_s^Lc_b^R}{c_s^Rc_b^L}\right)`", "`has_Nernst_diffusion_layer==True`" + "NDL thickness, cem", ":math:`\Delta^{L/R} = \frac{F D_b c_b^{L/R}}{\left(t_+^{iem}-t_+ \right) i_{lim}}`", "`has_Nernst_diffusion_layer==True`" + "NDL thickness, aem", ":math:`\Delta^{L/R} = - \frac{F D_b c_b^{L/R}}{\left(t_+^{iem}-t_+\right) i_{lim}}`", "`has_Nernst_diffusion_layer==True`" + "Concentration polarization ratio, cem", ":math:`\frac{c_s^L}{c_b^L} = 1+\frac{i}{i_{lim}},\qquad \frac{c_s^R}{c_b^R} = 1-\frac{i}{i_{lim}}`", "`has_Nernst_diffusion_layer==True` \ :sup:`1`" + "Concentration polarization ratio, aem", ":math:`\frac{c_s^L}{c_b^L} = 1-\frac{i}{i_{lim}},\qquad \frac{c_s^R}{c_b^R} = 1+\frac{i}{i_{lim}}`", "`has_Nernst_diffusion_layer==True`" **Note** diff --git a/docs/technical_reference/unit_models/electrodialysis_1D.rst b/docs/technical_reference/unit_models/electrodialysis_1D.rst index dcf320f061..33961a5f8e 100644 --- a/docs/technical_reference/unit_models/electrodialysis_1D.rst +++ b/docs/technical_reference/unit_models/electrodialysis_1D.rst @@ -257,7 +257,7 @@ This model can optionally calculate pressured drops along the flow path in the d "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`" + "Hydraulic diameter", ":math:`d_H=\frac{2db\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`" "Renold 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`" diff --git a/watertap/unit_models/electrodialysis_0D.py b/watertap/unit_models/electrodialysis_0D.py index ef0686f2ed..8da008024e 100644 --- a/watertap/unit_models/electrodialysis_0D.py +++ b/watertap/unit_models/electrodialysis_0D.py @@ -49,15 +49,15 @@ from watertap.core import ControlVolume0DBlock, InitializationMixin from watertap.costing.unit_models.electrodialysis import cost_electrodialysis -__author__ = " Xiangyu Bi, Austin Ladshaw," +__author__ = " Xiangyu Bi, Austin Ladshaw, Kejia Hu" _log = idaeslog.getLogger(__name__) class LimitingCurrentDensityMethod(Enum): InitialValue = 0 - # Empirical = 1 - # Theoretical = 2 TODO: 1 and 2 + Empirical = 1 + Theoretical = 2 class ElectricalOperationMode(Enum): @@ -1230,6 +1230,69 @@ def eq_current_dens_lim_ioa(self, t): for j in self.cation_set ) ) + elif ( + self.config.limiting_current_density_method + == LimitingCurrentDensityMethod.Theoretical + ): + self._get_fluid_dimensionless_quantities() + return self.current_dens_lim_ioa[t] == ( + self.N_Sh + * self.diffus_mass + * self.hydraulic_diameter**-1 + * Constants.faraday_constant + * ( + sum( + self.ion_trans_number_membrane["cem", j] + / self.config.property_package.charge_comp[j] + for j in self.cation_set + ) + - sum( + self.diluate.properties_in[t].trans_num_phase_comp["Liq", j] + / self.config.property_package.charge_comp[j] + for j in self.cation_set + ) + ) + ** -1 + * sum( + self.config.property_package.charge_comp[j] + * 0.5 + * ( + self.diluate.properties_in[t].conc_mol_phase_comp["Liq", j] + + self.diluate.properties_out[t].conc_mol_phase_comp[ + "Liq", j + ] + ) + for j in self.cation_set + ) + ) + elif ( + self.config.limiting_current_density_method + == LimitingCurrentDensityMethod.Empirical + ): + self.param_b = Param( + initialize=0.5, + units=pyunits.dimensionless, + doc="empirical parameter b to calculate limiting current density", + ) + self.param_a = Param( + initialize=25, + units=pyunits.coulomb + * pyunits.mol**-1 + * pyunits.meter ** (1 - self.param_b) + * pyunits.second ** (self.param_b - 1), + doc="empirical parameter a to calculate limiting current density", + ) + return self.current_dens_lim_ioa[ + t + ] == self.param_a * self.velocity_diluate[t] ** self.param_b * sum( + self.config.property_package.charge_comp[j] + * 0.5 + * ( + self.diluate.properties_in[t].conc_mol_phase_comp["Liq", j] + + self.diluate.properties_out[t].conc_mol_phase_comp["Liq", j] + ) + for j in self.cation_set + ) @self.Constraint( self.membrane_set, @@ -2143,15 +2206,6 @@ def calculate_scaling_factors(self): ), ) - if hasattr(self, "current_dens_lim_ioa"): - if iscale.get_scaling_factor(self.current_dens_lim_ioa) is None: - if ( - self.config.limiting_current_density_method - == LimitingCurrentDensityMethod.InitialValue - ): - sf = self.config.limiting_current_density_data**-1 - iscale.set_scaling_factor(self.current_dens_lim_ioa, sf) - if hasattr(self, "potential_nonohm_membrane_ioa"): if iscale.get_scaling_factor(self.potential_nonohm_membrane_ioa) is None: sf = ( @@ -2294,6 +2348,47 @@ def calculate_scaling_factors(self): ) iscale.set_scaling_factor(self.pressure_drop_total, sf) + if hasattr(self, "current_dens_lim_ioa"): + if iscale.get_scaling_factor(self.current_dens_lim_ioa) is None: + if ( + self.config.limiting_current_density_method + == LimitingCurrentDensityMethod.InitialValue + ): + sf = self.config.limiting_current_density_data**-1 + iscale.set_scaling_factor(self.current_dens_lim_ioa, sf) + elif ( + self.config.limiting_current_density_method + == LimitingCurrentDensityMethod.Empirical + ): + sf = 25**-1 * sum( + iscale.get_scaling_factor( + self.diluate.properties_in[0].conc_mol_phase_comp["Liq", j] + ) + ** 2 + for j in self.cation_set + ) ** (0.5 * 0.5) + iscale.set_scaling_factor(self.current_dens_lim_ioa, sf) + elif ( + self.config.limiting_current_density_method + == LimitingCurrentDensityMethod.Theoretical + ): + sf = ( + iscale.get_scaling_factor(self.N_Sh) + * iscale.get_scaling_factor(self.diffus_mass) + * iscale.get_scaling_factor(self.hydraulic_diameter) ** -1 + * 96485**-1 + * sum( + iscale.get_scaling_factor( + self.diluate.properties_in[0].conc_mol_phase_comp[ + "Liq", j + ] + ) + ** 2 + for j in self.cation_set + ) + ** 0.5 + ) + iscale.set_scaling_factor(self.current_dens_lim_ioa, sf) # Constraint scaling for ind, c in self.eq_current_voltage_relation.items(): iscale.constraint_scaling_transform( diff --git a/watertap/unit_models/tests/test_electrodialysis_0D.py b/watertap/unit_models/tests/test_electrodialysis_0D.py index 16b23ef69c..c0c94577cc 100644 --- a/watertap/unit_models/tests/test_electrodialysis_0D.py +++ b/watertap/unit_models/tests/test_electrodialysis_0D.py @@ -19,6 +19,7 @@ FrictionFactorMethod, HydraulicDiameterMethod, ) +from watertap.unit_models.electrodialysis_0D import LimitingCurrentDensityMethod from watertap.costing import WaterTAPCosting from pyomo.environ import ( ConcreteModel, @@ -1761,3 +1762,159 @@ def test_deltaP_configerr(self): pressure_drop_method=PressureDropMethod.Darcy_Weisbach, has_pressure_change=False, ) + + +class Test_Limiting_Current_Density_Method: + @pytest.fixture(scope="class") + def ed_l0(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-"], + "mw_data": {"H2O": 18e-3, "Na_+": 23e-3, "Cl_-": 35.5e-3}, + "elec_mobility_data": {("Liq", "Na_+"): 5.19e-8, ("Liq", "Cl_-"): 7.92e-8}, + "charge": {"Na_+": 1, "Cl_-": -1}, + "diffusivity_data": {("Liq", "Na_+"): 1.33e-9, ("Liq", "Cl_-"): 2.03e-9}, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Electrodialysis0D( + property_package=m.fs.properties, + operation_mode=ElectricalOperationMode.Constant_Voltage, + has_Nernst_diffusion_layer=True, + has_nonohmic_potential_membrane=True, + limiting_current_density_method=LimitingCurrentDensityMethod.InitialValue, + limiting_current_density_data=500, + ) + return m + + @pytest.fixture(scope="class") + def ed_l1(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-"], + "mw_data": {"H2O": 18e-3, "Na_+": 23e-3, "Cl_-": 35.5e-3}, + "elec_mobility_data": {("Liq", "Na_+"): 5.19e-8, ("Liq", "Cl_-"): 7.92e-8}, + "charge": {"Na_+": 1, "Cl_-": -1}, + "diffusivity_data": {("Liq", "Na_+"): 1.33e-9, ("Liq", "Cl_-"): 2.03e-9}, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Electrodialysis0D( + property_package=m.fs.properties, + operation_mode=ElectricalOperationMode.Constant_Voltage, + has_Nernst_diffusion_layer=True, + has_nonohmic_potential_membrane=True, + limiting_current_density_method=LimitingCurrentDensityMethod.Empirical, + ) + return m + + @pytest.fixture(scope="class") + def ed_l2(self): + m = ConcreteModel() + m.fs = FlowsheetBlock(dynamic=False) + ion_dict = { + "solute_list": ["Na_+", "Cl_-"], + "mw_data": {"H2O": 18e-3, "Na_+": 23e-3, "Cl_-": 35.5e-3}, + "elec_mobility_data": {("Liq", "Na_+"): 5.19e-8, ("Liq", "Cl_-"): 7.92e-8}, + "charge": {"Na_+": 1, "Cl_-": -1}, + "diffusivity_data": {("Liq", "Na_+"): 1.33e-9, ("Liq", "Cl_-"): 2.03e-9}, + } + m.fs.properties = MCASParameterBlock(**ion_dict) + m.fs.unit = Electrodialysis0D( + property_package=m.fs.properties, + operation_mode=ElectricalOperationMode.Constant_Voltage, + has_Nernst_diffusion_layer=True, + has_nonohmic_potential_membrane=True, + limiting_current_density_method=LimitingCurrentDensityMethod.Theoretical, + ) + return m + + @pytest.mark.unit + def test_limiting_current_various_methods(self, ed_l0, ed_l1, ed_l2): + ed_m = (ed_l0, ed_l1, ed_l2) + for m in ed_m: + m.fs.unit.water_trans_number_membrane["cem"].fix(5.8) + m.fs.unit.water_trans_number_membrane["aem"].fix(4.3) + m.fs.unit.water_permeability_membrane["cem"].fix(2.16e-14) + m.fs.unit.water_permeability_membrane["aem"].fix(1.75e-14) + m.fs.unit.voltage.fix(0.5) + m.fs.unit.electrodes_resistance.fix(0) + m.fs.unit.cell_pair_num.fix(10) + m.fs.unit.current_utilization.fix(1) + m.fs.unit.channel_height.fix(5e-4) + m.fs.unit.membrane_areal_resistance["cem"].fix(1.89e-4) + m.fs.unit.membrane_areal_resistance["aem"].fix(1.77e-4) + m.fs.unit.cell_width.fix(0.1) + m.fs.unit.cell_length.fix(0.79) + m.fs.unit.membrane_thickness["aem"].fix(1.3e-4) + m.fs.unit.membrane_thickness["cem"].fix(1.3e-4) + m.fs.unit.solute_diffusivity_membrane["cem", "Na_+"].fix(1.8e-10) + m.fs.unit.solute_diffusivity_membrane["aem", "Na_+"].fix(1.25e-10) + m.fs.unit.solute_diffusivity_membrane["cem", "Cl_-"].fix(1.8e-10) + m.fs.unit.solute_diffusivity_membrane["aem", "Cl_-"].fix(1.25e-10) + m.fs.unit.ion_trans_number_membrane["cem", "Na_+"].fix(1) + m.fs.unit.ion_trans_number_membrane["aem", "Na_+"].fix(0) + m.fs.unit.ion_trans_number_membrane["cem", "Cl_-"].fix(0) + m.fs.unit.ion_trans_number_membrane["aem", "Cl_-"].fix(1) + # set the inlet stream + m.fs.unit.inlet_diluate.pressure.fix(101325) + m.fs.unit.inlet_diluate.temperature.fix(298.15) + m.fs.unit.inlet_diluate.flow_mol_phase_comp[0, "Liq", "H2O"].fix(2.40e-1) + m.fs.unit.inlet_diluate.flow_mol_phase_comp[0, "Liq", "Na_+"].fix(7.38e-4) + m.fs.unit.inlet_diluate.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix(7.38e-4) + m.fs.unit.inlet_concentrate.pressure.fix(101325) + m.fs.unit.inlet_concentrate.temperature.fix(298.15) + m.fs.unit.inlet_concentrate.flow_mol_phase_comp[0, "Liq", "H2O"].fix( + 2.40e-1 + ) + m.fs.unit.inlet_concentrate.flow_mol_phase_comp[0, "Liq", "Na_+"].fix( + 7.38e-4 + ) + m.fs.unit.inlet_concentrate.flow_mol_phase_comp[0, "Liq", "Cl_-"].fix( + 7.38e-4 + ) + m.fs.unit.spacer_porosity.fix(1) + + m.fs.properties.set_default_scaling( + "flow_mol_phase_comp", 1e1, index=("Liq", "H2O") + ) + 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_-") + ) + + # Test ed_m0 + iscale.calculate_scaling_factors(ed_m[0]) + assert degrees_of_freedom(ed_m[0]) == 0 + initialization_tester(ed_m[0], outlvl=idaeslog.DEBUG) + badly_scaled_var_values = { + var.name: val for (var, val) in iscale.badly_scaled_var_generator(ed_m[0]) + } + assert not badly_scaled_var_values + results = solver.solve(ed_m[0]) + assert_optimal_termination(results) + assert value(ed_m[0].fs.unit.current_dens_lim_ioa[0]) == pytest.approx( + 444.002, + rel=1e-3, + ) + # Test ed_m1 + iscale.calculate_scaling_factors(ed_m[1]) + assert degrees_of_freedom(ed_m[1]) == 0 + initialization_tester(ed_m[1], outlvl=idaeslog.DEBUG) + results = solver.solve(ed_m[1]) + assert_optimal_termination(results) + assert value(ed_m[1].fs.unit.current_dens_lim_ioa[0]) == pytest.approx( + 353.041, rel=1e-3 + ) + # Test ed_m2 + iscale.calculate_scaling_factors(ed_m[2]) + ed_m[2].fs.unit.diffus_mass.fix(1.6e-9) + assert degrees_of_freedom(ed_m[2]) == 0 + initialization_tester(ed_m[2], outlvl=idaeslog.DEBUG) + results = solver.solve(ed_m[2]) + assert_optimal_termination(results) + assert value(ed_m[2].fs.unit.current_dens_lim_ioa[0]) == pytest.approx( + 281.100, rel=1e-3 + )