How to scale chemical reactions

Note

The following guide and examples assume that your GenericParameterBlock is named thermo_params, your GenericReactionParameterBlock is named rxn_params, your unit model is named unit, and your FlowsheetBlock is named fs. For information on setup of the basic chemistry modules, see How to setup simple chemistry.

No model in WaterTAP can solve without proper scaling of the constraints and variables within that model. This can be a difficult task for aqueous chemistry systems as these systems must deal with extremely dilute or trace chemical species and reaction coefficients that may vary from each other by several orders of magnitude.

To scale a set of chemical reactions properly requires different function arguments to be invoked depending on: (i) where the reactions are located in the model and (ii) what types of reactions are called for. Each of these different scenarios will be discussed and demonstrated below.

Note

IMPORTANT: Scaling of just the chemical reactions is insufficient for solving a chemistry module. User’s MUST also scale the chemical species in the system due to the dilute nature of aqueous chemistry. See How to scale chemical species. Additionally, user’s MUST also scale the energy balance equations. See also How to scale energy balance for chemistry.

After you have set all scaling factors for reactions, species, energy balances, etc., you need to call the ‘calculate_scaling_factors’ function on the entire model before attempting to solve. See below:

# After all scaling factors and constraint transformations are complete,
# call the following function on the model.
iscale.calculate_scaling_factors(model.fs.unit)

Types of Reactions

Inherent Reactions using log_power_law_equil form
Equilibrium Reactions using log_power_law_equil form
Rate Reactions using power_law_rate form
Solubility/Precipitation Reaction using log_solubility_product form
Stoichiometric Reactions

Inherent Reactions using log_power_law_equil form

Inherent reactions (see How to use inherent reactions) are placed into the GenericParameterBlock rather than the GenericReactionParameterBlock. As such, the scaling methods employed must be based on the thermo_config dictionary for this set of reactions. The code segment below shows a demonstration of scaling factors that generally work well for these types of reactions.

Inherent Scaling Demonstration

# Specify a minimum allowable scaling division factor
min_scale = 1e-3
#Add scaling factors for reactions by looping through the reactions
for i in model.fs.unit.control_volume.inherent_reaction_extent_index:
    # i[0] = time, i[1] = reaction

    # Grab the 'k_eq_ref' value directly from the thermo config dictionary
    scale = max(min_scale, thermo_config["inherent_reactions"][i[1]]["parameter_data"]["k_eq_ref"][0])
    iscale.set_scaling_factor(model.fs.unit.control_volume.inherent_reaction_extent[0.0,i[1]], 10/scale)
    iscale.constraint_scaling_transform(model.fs.unit.control_volume.properties_out[0.0].
            inherent_equilibrium_constraint[i[1]], 0.1)

What this code segment is doing is to first create scaling factors for the inherent_reaction_extent variable based on the inverse of the k_eq_ref value is (or based on some minimum inverse value).

Lastly, we perform a scaling transformation of the inherent_equilibrium_constraint. In our case, we scale this constraint down by a factor of 0.1 because the log_power_law_equil constraint form is already scaled up reasonably well.

Equilibrium Reactions using log_power_law_equil form

Equilibrium reactions (see How to setup simple chemistry) are placed into the GenericReactionParameterBlock. As such, the scaling methods employed must be applied to the rxn_config dictionary for this set of reactions. The code segment below shows a demonstration of scaling factors that generally work well for these types of reactions.

Equilibrium Scaling Demonstration

# Specify a minimum allowable scaling division factor
min_scale = 1e-3
#Add scaling factors for reactions
for i in model.fs.unit.control_volume.equilibrium_reaction_extent_index:
    # i[0] = time, i[1] = reaction

    # Grab the 'k_eq_ref' value from the reaction config
    scale = max(min_scale, rxn_config["equilibrium_reactions"][i[1]]["parameter_data"]["k_eq_ref"][0])
    iscale.set_scaling_factor(model.fs.unit.control_volume.equilibrium_reaction_extent[0.0,i[1]], 10/scale)
    iscale.constraint_scaling_transform(
        model.fs.unit.control_volume.reactions[0.0].equilibrium_constraint[i[1]], 0.1)

Note

These scaling arguments are identical to the Inherent Reaction scaling methods, however, because these reactions exist in a different location of the model, we showed this here for completeness. All reactions, regardless of location, need scaling.

Rate Reactions using power_law_rate form

Rate reactions only exist in the GenericReactionParameterBlock and so these scaling arguments apply to rxn_params for these types of reactions. These are much simpler to scale than both the Inherent and Equilibrium reactions, but are just as important to apply scaling for. Below is a demonstration of applying scaling.

Rate Reaction Scaling Demonstration

# Scaling for kinetic reactions
for i in model.fs.rxn_params.rate_reaction_idx:
    scale = value(model.fs.unit.control_volume.reactions[0.0].reaction_rate[i].expr)
    iscale.set_scaling_factor(model.fs.unit.control_volume.rate_reaction_extent[0.0,i], 1000/scale)

Note

We only need to call ‘set_scaling_factor’ here and NOT call ‘constraint_scaling_transform’ because this scaling factor will cascade into the constraints automatically once you call ‘calculate_scaling_factors’ on the model. This is different from other reactions because there is no ‘log form’ for rate reactions. The ‘log form’ always requires some additional treatment.

Solubility/Precipitation Reaction using log_solubility_product form

To scale these reactions, you will use the same methods outlined above for Equilibrium and Inherent reactions. However, there is an additional step. That additional step involves setting a smoothing parameter eps (which is a factor unique to the ‘log_solubility_product’ function). Below is a demonstration of setting up that smoothing parameter assuming your solubility reactions are in the rxn_params object and the rxn_config dictionary.

Setting eps Smoothing Factor for Solubility Products

# Specify a minimum allowable scaling factor for the eps
factor = 1e-2
for rid in model.fs.rxn_params.equilibrium_reaction_idx:
    # Grab the 'k_eq_ref' value from the reaction config
    scale = rxn_config["equilibrium_reactions"][rid]["parameter_data"]["k_eq_ref"][0]

    # NOTE: ONLY certain functions have an eps value that we need to set
    if (hasattr(rxn_params.component("reaction_"+rid),"eps")):
        # highest allowable value for setting eps based on k_eq_ref
        if scale < 1e-16:
            model.fs.rxn_params.component("reaction_"+rid).eps.value = scale*factor
        else:
            model.fs.rxn_params.component("reaction_"+rid).eps.value = 1e-16*factor

Stoichiometric Reactions

Stoichiometric reactions are generally the simplest to scale. However, determining how much to scale them by is not always clear. It depends on what are the expected changes in molar flows due to the reaction. Since there is not always a clear way to determine this, the demonstration below simply shows you where the scaling is applied to within the framework.

Setting Scaling Factor for Stoichiometric Reaction

The sample below just gives a brief demo of how to scale a stoichiometric reaction named ‘R1’ by a given factor. All stoichiometric reactions are located in the ‘control_volume’ of the model and named ‘rate_reaction_extent’. This is because Rate Reactions and Stoichiometric Reactions have very similar implementations in the IDAES framework.

# Specify a factor to scale by
factor = 1
iscale.set_scaling_factor(m.fs.unit.control_volume.rate_reaction_extent[0.0,'R1'], factor)