search

watertap-org / watertap

The WaterTAP development repository
https://watertap.readthedocs.io/en/latest
Other
59 stars 57 forks source link

High TDS feed brines in NF DSPM-DE block #1221

Closed AARON683 closed 9 months ago

AARON683 commented 9 months ago

I am investigating 3 different brines in the NF DSPM-DE block. They have the following compositions:

Brine 1: TDS - 60,500 ppm

H2O: 897.669 mol/s Na2+: 14.050 mol/s K+: 0.291 mol/s Mg_2+: 1.463 mol/s Ca2+: 0.208 mol/s Cl-: 16.659 mol/s SO4_2-: 0.843 mol/s

Brine 2: TDS - 140,000 ppm

H2O: 161.410 mol/s Na2+: 6.378 mol/s K+: 0.132 mol/s Mg_2+: 0.664 mol/s Ca2+: 0.094 mol/s Cl-: 7.261 mol/s SO4_2-: 0.383 mol/s

Brine 3: TDS - 220,000 ppm

H2O: 83.427 mol/s Na2+: 5.711 mol/s K+: 0.118 mol/s Mg_2+: 0.595 mol/s Ca2+: 0.084 mol/s Cl-: 6.503 mol/s SO4_2-: 0.343 mol/s

The model works fine when Brine 1 is used as an input and it produces good results. However when I try to run the model with the higher TDS brines (Brines 2 and 3) I get an InitializationError from the MCAS property package stating 'fs.unit.feed_side.properties_interface failed to initialize successfully. Please check the output logs for more information.' The output logs do not really provide any further information other than the NF block failed to initialize. I have a few thought on what might be potentially causing this error:

I would greatly appreciate any additional information regarding whether the issue stems from the mentioned possibilities or if it originates from a completely different cause.

adam-a-a commented 9 months ago

@AARON683 Thank you for raising this issue. Can you share your script or all model specifications? Also sharing the full traceback of the error would be helpful too. This issue could be related to certain flowsheet specifications (e.g., fixing variables to certain values, too stringent of a hardness constraint, etc.), inadequate scaling/initialization, or something else, but it's difficult to determine without knowing more details.

AARON683 commented 9 months ago

@adam-a-a Thank you for the response.

I use the following script for the basic model:

# Import Pyomo Packages
from pyomo.environ import (ConcreteModel, TransformationFactory, Objective, Var, Constraint, NonNegativeReals, 
    assert_optimal_termination, units as pyunits)
from pyomo.network import Arc
from pyomo.util.check_units import assert_units_consistent
from pyomo.environ import value

# Import IDAES Packages
from idaes.core import FlowsheetBlock
from idaes.core.solvers import get_solver
from idaes.core.util.initialization import propagate_state
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.models.unit_models import (Feed, Product, StateJunction)
import idaes.core.util.scaling as iscale

# Import WaterTap Packages
from watertap.unit_models.nanofiltration_DSPMDE_0D import (NanofiltrationDSPMDE0D, 
                                                           MassTransferCoefficient,
                                                           ConcentrationPolarizationType)
from watertap.property_models.multicomp_aq_sol_prop_pack import (MCASParameterBlock, 
                                                                 ActivityCoefficientModel,
                                                                 DensityCalculation)

# Import math
import math

# Create model
m = ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False)

# Parameters for property package taken from Labban et al.
m.fs.properties = MCASParameterBlock(
    solute_list=["Ca_2+",
                "Mg_2+",
                "SO4_2-", 
                "Na_+", 
                "Cl_-",
                "K_+",
                ],
    diffusivity_data={
        ("Liq", "Ca_2+"): 0.792e-9,
        ("Liq", "Mg_2+"): 0.706e-09,
        ("Liq", "SO4_2-"): 1.06e-9,
        ("Liq", "Na_+"): 1.33e-09,
        ("Liq", "Cl_-"): 2.03e-09,
        ("Liq", "K_+"): 1.957e-09,
    },
    mw_data={
        "H2O": 0.018,
        "Ca_2+":0.040,
        "Mg_2+":0.024,
        "SO4_2-": 0.096,
        "Na_+": 0.023,
        "Cl_-": 0.035,
        "K_+": 0.0391,
    },
    stokes_radius_data={
        "Ca_2+":3.09e-10,
        "Mg_2+":3.47e-10,
        "SO4_2-": 2.30e-10,
        "Cl_-": 1.21e-10,
        "Na_+": 1.84e-10,
        "K_+": 1.25e-10,
    },
    charge={
        "Ca_2+": 2,
        "Mg_2+": 2,
        "SO4_2-": -2,
        "Na_+": 1,
        "Cl_-": -1,
        "K_+": 1,
    },
    activity_coefficient_model=ActivityCoefficientModel.davies,
    density_calculation=DensityCalculation.constant,
)

# Add NF unit
m.fs.unit = NanofiltrationDSPMDE0D(property_package=m.fs.properties, 
                                mass_transfer_coefficient=MassTransferCoefficient.spiral_wound,
                                concentration_polarization_type=ConcentrationPolarizationType.calculated
                                )

# Number of vessels in parallel
vessels_num = 1

# Specify molar flow in mol/s
# The values specified is the total molar flow of brine produced in SP40 design
# This is divided between a number of pressure vessels in parallel (vessels_num)
molar_flow = {
                        "H2O": 161.410/vessels_num,
                        "Ca_2+": 0.094/vessels_num,
                        "Mg_2+": 0.664/vessels_num,
                        "SO4_2-": 0.383/vessels_num,
                        "Cl_-": 7.562/vessels_num,
                        "Na_+": 6.378/vessels_num,
                        "K_+": 0.132/vessels_num,
                        }

# Fix the molar flowrates of each component
for ion, x in molar_flow.items():
    mol_comp_flow = x
    m.fs.unit.inlet.flow_mol_phase_comp[0, "Liq", ion].fix(mol_comp_flow)

# Scaling of feed flowrates
m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1e-2, index=("Liq", "H2O"))
m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1e2, index=("Liq", "Ca_2+"))
m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1e1, index=("Liq", "Mg_2+"))
m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1e1, index=("Liq", "SO4_2-"))
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", "Na_+"))
m.fs.properties.set_default_scaling("flow_mol_phase_comp", 1e1, index=("Liq", "K_+"))
iscale.set_scaling_factor(m.fs.unit.area, 1e-1)
iscale.calculate_scaling_factors(m.fs)

# Use assert electroneutrality method from property model to ensure the ion concentrations provided obey electroneutrality condition
m.fs.unit.feed_side.properties_in[0].assert_electroneutrality(
tee=False,
defined_state=True,
adjust_by_ion="Cl_-",
get_property="flow_mol_phase_comp",
)

# Specify inlet state variables
m.fs.unit.inlet.temperature[0].fix(273.15+37)               # DEWA Brine Temperature
m.fs.unit.inlet.pressure[0].fix((20*1e5) + 101325)          # TMP drop + atomospheric pressure = inlet pressure (Pa)

# Fix membrane properties
m.fs.unit.radius_pore.fix(0.507e-9)                         # Pore radius of membrane (m) Micari et al.
m.fs.unit.membrane_thickness_effective.fix(0.8e-6)          # Membrane thickness (m) Micari et al.
m.fs.unit.membrane_charge_density.fix(-50)                  # Charge density (mol/m3) Artificial Seawater in Micari et al.
m.fs.unit.dielectric_constant_pore.fix(42.5)                # Dielectric constant (-) Micari et al.

# Fix permeate pressure at atomospheric
m.fs.unit.permeate.pressure[0].fix(101325)                  # Permeate pressure (Pa)

# Fix additional membrane properties
m.fs.unit.spacer_porosity.fix(0.85)                         # Spacer porosity (-) 
m.fs.unit.channel_height.fix(7e-4)                          # Channel height (m) Roy et al.
m.fs.unit.velocity[0, 0].fix(0.3)                           # Cross flow velocity (m/s) Labban et al.
m.fs.unit.area.fix(7.6)                                     # Membrane area (m2) Dupont NF270 element
m.fs.unit.spacer_mixing_efficiency.fix(1)                   # Mixing efficiency (-) assumed completly mixed
m.fs.unit.spacer_mixing_length.fix(0.006)                   # Mixing length (m) Senthilmurugan et al. 
                                                            # Intialized as 6m in source code (think this is an error) 

# Intialize and solve the system
solver = get_solver()
m.fs.unit.initialize(optarg=solver.options)
assert_units_consistent(m)
assert degrees_of_freedom(m) == 0
solver.solve(m, tee=False)

# Return the results of model
m.fs.unit.report()

This model runs fine when I use Brine 1 as the inlet feed. However when I use Brine 2 and Brine 3, which have higher TDS. Here is the full traceback of they error:

---------------------------------------------------------------------------
InitializationError                       Traceback (most recent call last)
[~\AppData\Local\Temp\ipykernel_29212\2009364540.py](https://file+.vscode-resource.vscode-cdn.net/c%3A/Users/Aaron/OneDrive/Documents/Chemical%20Engineering%20YEAR%204/Reserach%20Project%20-%20Desolenator/Thesis/4-Nanofiltration/NF%20Python%20Modelling/~/AppData/Local/Temp/ipykernel_29212/2009364540.py) in <cell line: 151>()
    149 # Intialize and solve the system
    150 solver = get_solver()
--> 151 m.fs.unit.initialize(optarg=solver.options)
    152 assert_units_consistent(m)
    153 assert degrees_of_freedom(m) == 0

[c:\Users\Aaron\anaconda3\envs\Water\lib\site-packages\watertap\core\initialization_mixin.py](file:///C:/Users/Aaron/anaconda3/envs/Water/lib/site-packages/watertap/core/initialization_mixin.py) in initialize(self, *args, **kwargs)
     21     def initialize(self, *args, **kwargs):
     22         try:
---> 23             return super().initialize(*args, **kwargs)
     24         except InitializationError:
     25             for blk in self._initialization_order:

[c:\Users\Aaron\anaconda3\envs\Water\lib\site-packages\idaes\core\base\unit_model.py](file:///C:/Users/Aaron/anaconda3/envs/Water/lib/site-packages/idaes/core/base/unit_model.py) in initialize(blk, *args, **kwargs)
    528 
    529         # Remember to collect flags for fixed vars
--> 530         flags = blk.initialize_build(*args, **kwargs)
    531 
    532         # If costing block exists, activate and initialize

[c:\Users\Aaron\anaconda3\envs\Water\lib\site-packages\watertap\unit_models\nanofiltration_DSPMDE_0D.py](file:///C:/Users/Aaron/anaconda3/envs/Water/lib/site-packages/watertap/unit_models/nanofiltration_DSPMDE_0D.py) in initialize_build(self, initialize_guess, state_args, outlvl, solver, optarg, fail_on_warning, ignore_dof, automate_rescale)
   1387             state_args=state_args["retentate"],
   1388         )
-> 1389         self.feed_side.properties_interface.initialize(
   1390             outlvl=outlvl,
   1391             optarg=optarg,

[c:\Users\Aaron\anaconda3\envs\Water\lib\site-packages\watertap\property_models\multicomp_aq_sol_prop_pack.py](file:///C:/Users/Aaron/anaconda3/envs/Water/lib/site-packages/watertap/property_models/multicomp_aq_sol_prop_pack.py) in initialize(self, state_args, state_vars_fixed, hold_state, outlvl, solver, optarg)
    810 
    811         if (not skip_solve) and (not check_optimal_termination(results)):
--> 812             raise InitializationError(
    813                 f"{self.name} failed to initialize successfully. Please "
    814                 f"check the output logs for more information."

InitializationError: fs.unit.feed_side.properties_interface failed to initialize successfully. Please check the output logs for more information.
adam-a-a commented 9 months ago

@AARON683 Thank you again for raising this issue. Your probing is very useful, and we highly encourage you to continue doing so and posting issues as they arise. I noticed that you identified a subtle error in the intended initial value for spacer_mixing_length--it was meant to be 0.006m. I will update this although the user could arguably provide his/her own value for this. Now- let's move to the real issue at hand.

The reason why the properties_interface stateblock failed to initialize was because of the calculation of activity coefficients (act_coeff_phase_comp) and the fact that an upper bound of 1.001 is set as the default upper bound for act_coeff_phase_comp in the MCAS property model.

Setting the upper bound for act_coeff_phase_comp to None on that stateblock and subsequently running your script would then result in successfully initializing the properties_interface stateblock followed by a failure in initialization of the next stateblock that uses act_coeff_phase_comp. However, when setting the upper bound of act_coeff_phase_comp to None globally, the model solves. Overall, I tried this for your second and third brine compositions, and both cases solved when I eliminated bounds on act_coeff_phase_comp. Checking the results for act_coeff_phase_comp, the highest value was something around ~ 1.02, indicating that there isn't really a need for setting an upper bound in the first place. I plan to remove bounds on act_coeff_phase_comp in a subsequent PR.

Another important note to consider is that MCAS currently supports the use of ideal activity coefficients (act_coeff_phase_comp=1) or use of the Davies model to compute activity coefficients. Thus, another solution to your issue would be to set the following configuration option to activity_coefficient_model=ActivityCoefficientModel.ideal instead of activity_coefficient_model=ActivityCoefficientModel.davies. More importantly though, the Davies model is applicable for ionic strengths up to 0.5M, and I think the ionic strengths of your "brine 2" and "brine 3" compositions exceed the valid range for the Davies model. Additionally, bear in mind that there may be other complexities that could limit the validity of the DSPM-DE model for high-salinity applications (e.g., Foo et al., 2023).

In summary, the restrictive upper bound on the act_coeff_phase_comp variable caused initialization/solve failures and removing said bound resolves the issue. However, be wary of applying the model to high-salinity applications since (1) activity coefficient models that would be valid for high salinities are not yet supported (e.g., Pitzer) and (2) the DSPM-DE model itself has limitations under such conditions that one should keep in mind.

EDIT: I should mention that with regard to activity coefficients, WaterTAP does have the capability to use OLI to compute activity coefficients for high-salinity compositions, but one would need an OLI license to make use of that extended functionality. Furthermore, we are working on more seamlessly integrating property models (e.g., MCAS) with the OLI Cloud API so that one could obtain properties that might be more difficult to calculate for an arbitrary feed composition. Alternatively, there might be other tools under consideration that we can link with WaterTAP property models (e.g., PHREEQC) which could be used to calculate activity coefficients at high salinities.

AARON683 commented 9 months ago

@adam-a-a Thank you for the comprehensive response. I am glad to hear that the true value of the spacer_mixing_length is supposed to be 0.006 m. While I was checking the results produced by model I created using the NF DSPM-DE block replicated those in literature (Labban et al., Micari et al.), I noticed the trends weren't matching up. However, once the value of the spacer_mixing_length was change to 0.006 m, the trends were very similar with minimal errors', hence validating the model.

image Permeate flux vs Ion Rejection rate for Artificial seawater solution – Labban et al. image Permeate flux vs Ion Rejection rate for Artificial seawater solution - WaterTap Model

Regarding the error with high TDS brines, I did not consider the activity coefficient model as the source of error, so thank you for this! I shall change investigate using the 'ideal' model and also setting the upper bound to none. Thank you for pointing out that both thew Davies activity coefficient model and the DSPM-DE model itself to not provide accurate predictions for high salinity feed solutions - I will keep this in mind during further analysis.

I look forward to the new additions you are working on!

adam-a-a commented 9 months ago

@AARON683 Excellent work validating against the literature! This feedback is great, and we look forward to future discussions with you on WaterTAP issues!