Pythonize Regressed Physical Chemical Propeties in DWSIM

Originally posted on

Initiating roadmap #

Physical properties like viscosity, density, thermal conductivity play a major role in accessing fluid conveyance power requirements as well as thermal requirements in a process thus properly predicting these can move whole simulation near to reality.

There is a python tutorial for overriding physical properties of material streams in DWSIM but how do to get the mathematical function that relates this physical property to state properties (Temperature / Pressure) of material stream as well as its composition if it’s a mixture. This is especially important if you are modelling strong ionic species like caustics or brines.

In this tutorial I will turn a graph that displays aqueous NaOH viscosity as a function of temperature & solute mass%.

1. Acquire data or graph from a reliable source. Eyeball the graph to generate a data table
2. Use multivariate nonlinear regression to get an approximate mathematical function
3. Pythonize that mathematical function for DWSIM to use

You can (& you should) use the same steps for any other physical property. For step 3] you can use Microsoft Excel or WPS / LibreCalc (if you are on Linux). You can also use python numpy for nonlinear regression as explained here.

Acquiring Data #

NaOH viscosity as a function of temperature & solute mass%

Regressing Data #

First you need to plot a curve relating x2 variables while keeping the third constant. Use polynomial (order=1,2 or 3), log, ln relationships to visualize which relationship suits best (R^2 value close to 1).

Viscosity & Temperature at 30% NaOH w.t%

Clever use of INDEX() & LINEST() functions can dynamically generate mathematical constants (a & b). You can read more about it here.

Now we need to so second regression that is variation of a & b with respect to NaOH mass fractions because after all viscosity of brines are dependent on x2 variables (temperature & solute mass fraction)

This completes our analysis & we are able to mathematically model viscosity using x3 math functions $$ \nu = \left| \begin{array}{l} f_1 (a, b, T) \\ f_2 (a, m_F) \\ f_3 (b, m_F) \end{array} \right. $$ We can now generate a table using above mathematical model & compare with experimental data.

As one can see regressed mathematical model is properly predicting viscosity especially at higher solute mass fractions

Pythonise function #

We can now pythonize this model in DWSIM & it will overwrite calculated viscosity which was highly underspecified.

# Multiple regression functions
import math as m
def calA(m_fr):
        return 5.7357 * m.exp(17.717*m_fr)

def calB(m_fr):
        return -3.8345*m_fr**2 - 0.8307*m_fr - 0.7612

# properties are nullable values, so we need to check if they have a value first.
viscobj = matstr.GetPhase('OverallLiquid').Properties.viscosity # in Pa.s (cP / 1000)
if (viscobj <> None):
    currval = viscobj
    m_fr = matstr.GetPhase('OverallLiquid').Compounds['Sodium Hydroxide'].MassFraction
    # check if NaOH is present
    if (m_fr > 0):
        flowsheet.WriteMessage('current liquid viscosity (cp): ' + str(currval*1000))
        m_temp = matstr.GetTemperature() - 273.15 # in degC
        a = calA(m_fr) # calculate constant a  
        b = calB(m_fr) # calculate constant b
        propval = a * m_temp**b / 1000
        flowsheet.WriteMessage('calculated liquid viscosity (cp): ' + str(propval*1000))
    else:
        propval = currval
else:
    propval = 0.0 

[11/26/2023 9:52:51 PM] The flowsheet was calculated successfully. [0.1266s] [11/26/2023 9:52:54 PM] The flowsheet is being calculated, please wait… [11/26/2023 9:52:54 PM] The flowsheet was calculated successfully. [0.03743s] [11/26/2023 9:55:06 PM] current liquid viscosity (cp): 0.547113119959 [11/26/2023 9:55:06 PM] calculated liquid viscosity (cp): 7.85757761474

As you can see the default viscosity was 0.547 & the proper one is 7.85 & comparing with experimental data it is very near.

Further findings #

While following above steps I also modified density calculations but there is a catch. For example, I am inputting volumetric flowrate for NaOH (aqueous) thus mass flow is calculated using default density which is grossly underreported (either using Rackett or Costald). So once you have overwritten density you also have to recalculate mass flow ONLY for feed stream. Please read this forum where I appreciated Daniel Medeiros help:

Thread for density & mass flow update through scripting