diff --git a/docs/technical_reference/unit_models/stoichiometric_reactor.rst b/docs/technical_reference/unit_models/stoichiometric_reactor.rst index 26b1f0bf39..2c3ca12cd7 100644 --- a/docs/technical_reference/unit_models/stoichiometric_reactor.rst +++ b/docs/technical_reference/unit_models/stoichiometric_reactor.rst @@ -133,9 +133,11 @@ Variables "Reagent dose", reagent_dose,[reagent],kg/:math:`\text{m}^3` "Reagent density", density_reagent,[reagent],kg/:math:`\text{m}^3` "Reagent flow mass", flow_mass_reagent,[reagent],kg/s + "Reagent molar flow", flow_mol_reagent,[reagent],mol/s "Reagent flow volume", flow_vol_reagent,[reagent],:math:`\text{m}^3`/s "Stoichiometric coefficients for dissolution", dissolution_stoich_comp, "[reagent, :math:`j`]",dimensionless "Flow mass of precipitant",flow_mass_precipitate,[precipitant],kg/s + "Molar flow of precipitant",flow_mol_precipitate,[precipitant],mol/s "Mass concentration of precipitant",conc_mass_precipitate,[precipitant],kg/:math:`\text{m}^3` "Stoichiometric coefficients for precipitation", precipitation_stoich_comp, "[precipitant, :math:`j`]",dimensionless "Fraction of solids in waste stream", waste_mass_frac_precipitate, None, fraction diff --git a/watertap/unit_models/stoichiometric_reactor.py b/watertap/unit_models/stoichiometric_reactor.py index 262c576a53..dedbb1ec50 100644 --- a/watertap/unit_models/stoichiometric_reactor.py +++ b/watertap/unit_models/stoichiometric_reactor.py @@ -16,7 +16,6 @@ Var, Param, Suffix, - NonNegativeReals, Reals, Reference, TransformationFactory, @@ -281,30 +280,28 @@ def build(self): ) # Add unit parameters if self.has_dissolution_reaction: - self.mw_reagent = Var( + self.mw_reagent = Param( self.reagent_list, units=pyunits.kg / pyunits.mol, doc="Molecular weight of reagents", + mutable=True, ) - self.density_reagent = Var( + self.density_reagent = Param( self.reagent_list, initialize=1000, units=units_meta("mass") / units_meta("volume"), doc="Density of reagents", + mutable=True, ) for r in self.reagent_list: - self.mw_reagent[r].fix( - pyunits.convert( - self.config.reagent[r]["mw"], - to_units=pyunits.kg / pyunits.mol, - ) + self.mw_reagent[r] = pyunits.convert( + self.config.reagent[r]["mw"], + to_units=pyunits.kg / pyunits.mol, ) if self.config.reagent[r].get("density_reagent") is not None: - self.density_reagent[r].fix( - pyunits.convert( - self.config.reagent[r].get("density_reagent"), - to_units=(units_meta("mass") / units_meta("volume")), - ) + self.density_reagent[r] = pyunits.convert( + self.config.reagent[r].get("density_reagent"), + to_units=(units_meta("mass") / units_meta("volume")), ) else: _log.warning( @@ -312,11 +309,9 @@ def build(self): r ) ) - self.density_reagent[r].fix( - pyunits.convert( - 1 * pyunits.kg / pyunits.liter, - to_units=(units_meta("mass") / units_meta("volume")), - ) + self.density_reagent[r] = pyunits.convert( + 1 * pyunits.kg / pyunits.liter, + to_units=(units_meta("mass") / units_meta("volume")), ) self.dissolution_stoich_comp = Param( self.reagent_list, @@ -334,21 +329,28 @@ def build(self): self.reagent_dose = Var( self.reagent_list, initialize=1, - domain=NonNegativeReals, + domain=Reals, units=units_meta("mass") / units_meta("volume"), doc="reagent dose", ) self.flow_mass_reagent = Var( self.reagent_list, initialize=1e-3, - domain=NonNegativeReals, + domain=Reals, units=units_meta("mass") / units_meta("time"), doc="Mass flowrate of reagent", ) + self.flow_mol_reagent = Var( + self.reagent_list, + initialize=1e-3, + domain=Reals, + units=pyunits.mol / units_meta("time"), + doc="Mass flowrate of reagent", + ) self.flow_vol_reagent = Var( self.reagent_list, initialize=1e-3, - domain=NonNegativeReals, + domain=Reals, units=units_meta("volume") / units_meta("time"), doc="Volume flowrate of reagent", ) @@ -368,17 +370,16 @@ def build(self): units=pyunits.dimensionless, doc="Solid mass fraction in sludge", ) - self.mw_precipitate = Var( + self.mw_precipitate = Param( self.precipitate_list, units=pyunits.kg / pyunits.mol, doc="Molecular weight of precipitate", + mutable=True, ) for p in self.precipitate_list: - self.mw_precipitate[p].fix( - pyunits.convert( - self.config.precipitate[p]["mw"], - to_units=pyunits.kg / pyunits.mol, - ) + self.mw_precipitate[p] = pyunits.convert( + self.config.precipitate[p]["mw"], + to_units=pyunits.kg / pyunits.mol, ) self.precipitation_stoich_comp = Param( @@ -399,14 +400,21 @@ def build(self): self.flow_mass_precipitate = Var( self.precipitate_list, initialize=1e-3, - domain=NonNegativeReals, + domain=Reals, units=units_meta("mass") / units_meta("time"), doc="Mass flowrate of precipitate", ) + self.flow_mol_precipitate = Var( + self.precipitate_list, + initialize=1e-3, + domain=Reals, + units=pyunits.mol / units_meta("time"), + doc="Mass flowrate of precipitate", + ) self.conc_mass_precipitate = Var( self.precipitate_list, initialize=1, - domain=NonNegativeReals, + domain=Reals, units=units_meta("mass") / units_meta("volume"), doc="Mass concentration of precipitate", ) @@ -585,6 +593,14 @@ def eq_dissolution_reaction_generation(b, t, j): for r in self.reagent_list ) + @self.Constraint( + self.flowsheet().config.time, + self.reagent_list, + doc="Mole generation from dissolution", + ) + def eq_flow_mol_reagent(b, t, r): + return b.flow_mass_reagent[r] / b.mw_reagent[r] == b.flow_mol_reagent[r] + @self.Constraint( self.flowsheet().config.time, self.reagent_list, @@ -604,8 +620,8 @@ def eq_flow_mass_reagent(b, t, r): ) def eq_flow_vol_reagent(b, t, r): return ( - b.flow_mass_reagent[r] - == b.flow_vol_reagent[r] * b.density_reagent[r] + b.flow_mass_reagent[r] / b.density_reagent[r] + == b.flow_vol_reagent[r] ) if self.has_precipitation_reaction: @@ -649,6 +665,17 @@ def eq_precipitation_reaction_generation(b, t, j): for p in self.precipitate_list ) + @self.Constraint( + self.flowsheet().config.time, + self.precipitate_list, + doc="Mole generation from precipitation", + ) + def eq_flow_mol_precipitate(b, t, p): + return ( + b.flow_mass_precipitate[p] / b.mw_precipitate[p] + == b.flow_mol_precipitate[p] + ) + @self.Constraint( self.flowsheet().config.time, self.precipitate_list, @@ -784,7 +811,10 @@ def calculate_scaling_factors(self): # these variables should have user input, if not there will be a warning if self.has_dissolution_reaction: + for r in self.reagent_list: + sf = 1 / self.mw_reagent[r].value + iscale.set_scaling_factor(self.mw_reagent[r], sf) if iscale.get_scaling_factor(self.reagent_dose[r]) is None: if iscale.get_scaling_factor(self.flow_mass_reagent[r]) is not None: # scale reagent_dose based on flow_mass_reagent @@ -802,15 +832,10 @@ def calculate_scaling_factors(self): self.reagent_dose[r], default=1, warning=True ) iscale.set_scaling_factor(self.reagent_dose[r], sf) - if iscale.get_scaling_factor(self.flow_mass_reagent[r]) is None: + if iscale.get_scaling_factor(self.flow_vol_reagent[r]) is None: sf = iscale.get_scaling_factor( - self.dissolution_reactor.properties_in[0].flow_vol_phase["Liq"] + self.flow_mass_reagent[r], default=1e4, warning=True ) - sf = iscale.get_scaling_factor(self.reagent_dose[r]) * sf - iscale.set_scaling_factor(self.flow_mass_reagent[r], sf) - if iscale.get_scaling_factor(self.flow_vol_reagent[r]) is None: - sf = iscale.get_scaling_factor(self.flow_mass_reagent[r]) - print(sf, self.density_reagent[r].value) sf = sf / self.density_reagent[r].value iscale.set_scaling_factor(self.flow_vol_reagent[r], sf) for (t, j), v in self.dissolution_reaction_generation_comp.items(): @@ -855,9 +880,23 @@ def calculate_scaling_factors(self): self.dissolution_reactor.properties_in[0].flow_vol_phase["Liq"] ) iscale.set_scaling_factor(self.flow_mass_reagent[r], sf) + if iscale.get_scaling_factor(self.flow_mol_reagent[r]) is None: + sf = iscale.get_scaling_factor(self.flow_mass_reagent[r]) + sf = sf / iscale.get_scaling_factor(self.mw_reagent[r]) + iscale.set_scaling_factor(self.flow_mol_reagent[r], sf) + for (t, r), con in self.eq_flow_mol_reagent.items(): + sf = iscale.get_scaling_factor( + self.flow_mol_reagent[r], default=1, warning=True + ) + iscale.constraint_scaling_transform(con, sf) + iscale.constraint_scaling_transform( + self.dissolution_reactor.eq_isothermal_dissolution[0], 1e-2 + ) if self.has_precipitation_reaction: for p in self.precipitate_list: + sf = 1 / self.mw_precipitate[p].value + iscale.set_scaling_factor(self.mw_precipitate[p], sf) if iscale.get_scaling_factor(self.flow_mass_precipitate[p]) is None: if ( iscale.get_scaling_factor(self.conc_mass_precipitate[p]) @@ -878,7 +917,21 @@ def calculate_scaling_factors(self): self.flow_mass_precipitate[p], default=1e3, warning=True ) iscale.set_scaling_factor(self.flow_mass_precipitate[p], sf) + if iscale.get_scaling_factor(self.flow_mass_precipitate[p]) is None: + iscale.set_scaling_factor(self.flow_mass_precipitate[p], sf) + if iscale.get_scaling_factor(self.flow_mol_precipitate[p]) is None: + sf = iscale.get_scaling_factor(self.flow_mass_precipitate[p]) + sf = sf / iscale.get_scaling_factor(self.mw_precipitate[p]) + iscale.set_scaling_factor(self.flow_mol_precipitate[p], sf) + iscale.constraint_scaling_transform( + self.precipitation_reactor.eq_isothermal_precipitation[0], 1e-2 + ) + for (t, p), con in self.eq_flow_mol_precipitate.items(): + sf = iscale.get_scaling_factor( + self.flow_mol_precipitate[p], default=1, warning=True + ) + iscale.constraint_scaling_transform(con, sf) for p in self.precipitate_list: if iscale.get_scaling_factor(self.conc_mass_precipitate[p]) is None: sf = iscale.get_scaling_factor( diff --git a/watertap/unit_models/tests/test_stoichiometric_reactor.py b/watertap/unit_models/tests/test_stoichiometric_reactor.py index 1cc98822ca..46c68cc9f0 100644 --- a/watertap/unit_models/tests/test_stoichiometric_reactor.py +++ b/watertap/unit_models/tests/test_stoichiometric_reactor.py @@ -462,6 +462,10 @@ def test_molar(self, basic_unit_molar): m.fs.unit.flow_mass_reagent["CaO"] / (56.0774 * 1e-3) ) + assert pytest.approx(flow_ca_in_reagent, rel=1e-5) == value( + m.fs.unit.flow_mol_reagent["CaO"] + ) + assert pytest.approx(expected_ca_mol_flow, rel=1e-5) == value( m.fs.unit.dissolution_reactor.properties_out[0].flow_mol_phase_comp[ "Liq", "Ca_2+" @@ -617,7 +621,9 @@ def test_precipitation(self, precipitation_reactor): "Liq", "Ca_2+" ] ) - + assert pytest.approx(1e-3 / (100.09 * 1e-3), rel=1e-5) == value( + m.fs.unit.flow_mol_precipitate["Calcite"] + ) expected_mg_mol_flow = 0.1 - 1e-4 / (58.3197 * 1e-3) assert pytest.approx(expected_mg_mol_flow, rel=1e-5) == value(