Skip to content

Commit

Permalink
Add limiting current density methods to electrodialysis_0D (watertap-…
Browse files Browse the repository at this point in the history 
…org#1242)

* add a sim file and modified ED0D

* modified ED and sim files

* add files for testing 0D ED pressure drops

* validate

* modified test files for 1d and 0d

* revised on test 1d and 0d files

* delete extra files

* remove extra files

* clean the files

* modified the test values

* none

* blank the files

* 1

* debug the pylint

* modified the limiting current density

* add a practice file

* added tests for limiting current density methods

* remove files

* modify the limiting current density docs

* modify docs

* formatting

* modify equations

* resolve conflicts

* resolving

* resolving

* resolving

* resolve

* resolve

* resolve

* eq restructuring on one constr

* Remove unnecessary terminal print

---------

Co-authored-by: Xiangyu B <97477152+lbibl@users.noreply.github.com>
Co-authored-by: lbibl <xiangyu.bi@gmail.com>
Co-authored-by: Hunter Barber <101219154+hunterbarber@users.noreply.github.com>
Co-authored-by: Keith Beattie <ksbeattie@lbl.gov>
  • Loading branch information
@kejia6002 @lbibl
@hunterbarber @ksbeattie
5 people authored Jan 4, 2024
1 parent 65f9749 commit 8cff71f
Show file tree
Hide file tree
Showing 4 changed files with 275 additions and 20 deletions.
17 changes: 10 additions & 7 deletions docs/technical_reference/unit_models/electrodialysis_0D.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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**
Expand Down
2 changes: 1 addition & 1 deletion docs/technical_reference/unit_models/electrodialysis_1D.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`"
Expand Down
119 changes: 107 additions & 12 deletions watertap/unit_models/electrodialysis_0D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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 = (
Expand Down Expand Up @@ -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(
Expand Down
157 changes: 157 additions & 0 deletions watertap/unit_models/tests/test_electrodialysis_0D.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
FrictionFactorMethod,
HydraulicDiameterMethod,
)
from watertap.unit_models.electrodialysis_0D import LimitingCurrentDensityMethod
from watertap.costing import WaterTAPCosting
from pyomo.environ import (
ConcreteModel,
Expand Down Expand Up @@ -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
)

0 comments on commit 8cff71f

Please sign in to comment.