Fit Pure Fluid Parameters¶

This example shows how to use the new fitting class of version 0.21 of teqp to fit the model parameters for the SAFT-VR-Mie model based on fitting pseudo-experimental data obtained from a reference equation of state.

In version 0.23, the parameter optimization objects were moved into their own submodule called paramopt

[1]:

import teqp
import numpy as np
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP
import scipy.optimize
display(teqp.__version__)

nonpolar = {
    "kind": "SAFT-VR-Mie",
    "model": {
        "coeffs": [
        {
          "name": "R32",
          "BibTeXKey": "Bell",
          "m": 1.2476268271391935,
          "sigma_m": 3.6080717234117107e-10,
          "epsilon_over_k": 172.53065054286867,
          "lambda_r": 14.634722358167384,
          "lambda_a": 6
        }
      ]
    }
}
template = {
  'kind': 'genericSAFT',
  'model': {
      'nonpolar': nonpolar
  }
}

# NOTE: '/' inside field names MUST be escaped as ~1; see https://datatracker.ietf.org/doc/html/rfc6901#section-3
pointers = [
    '/model/nonpolar/model/coeffs/0/m',
    '/model/nonpolar/model/coeffs/0/sigma_m',
    '/model/nonpolar/model/coeffs/0/epsilon_over_k',
    '/model/nonpolar/model/coeffs/0/lambda_r',
    '/model/nonpolar/model/coeffs/0/lambda_a'
]
x0 = [1.5, 3e-10, 150, 19, 5.7]
bounds = [(1,5), (2e-10,5e-10), (100,400), (12,50), (5.1, 6.0)]

'0.22.0'

[2]:

FLD = 'Methane'
ppo = teqp.paramopt.PureParameterOptimizer(template, pointers)

Ts = np.linspace(100, 170, 10)

# Generate some artificial pseudo-experimental data from
# the reference EOS as implemented in CoolProp
for T in Ts:
    pt = teqp.paramopt.SatRhoLPWPoint()
    rhoL, rhoV = [CP.PropsSI('Dmolar','Q',Q,'T',T,FLD) for Q in [0,1]]
    p = CP.PropsSI('P','Q',0,'T',T,FLD)
    w = CP.PropsSI('speed_of_sound','Q',0,'T',T,FLD)
    pt.T = T
    # Measurands (here, pseudo-experimental values)
    pt.p_exp = p
    pt.rhoL_exp = rhoL
    pt.w_exp = w

    # Additional parameters
    pt.rhoL_guess = rhoL
    pt.rhoV_guess = rhoV
    pt.R = 8.31446261815324
    pt.M = CP.PropsSI('molemass',FLD)
    AS = CP.AbstractState('HEOS',FLD)
    AS.update(CP.DmolarT_INPUTS, 1e-10, T)
    pt.Ao20 = AS.tau()**2*AS.d2alpha0_dTau2() # -cv0/R

    # Weights (multiplied by 100 to put on a percentage basis)
    pt.weight_p = 1*100
    pt.weight_rho = 1*100
    pt.weight_w = 0.25*100
    ppo.add_one_contribution(pt)

def cost_function(x):
#     return ppo.cost_function_threaded(x, 10) # This is an option if you have lots of threads
    return ppo.cost_function(x)

r = scipy.optimize.differential_evolution(cost_function, bounds=bounds, disp=True, maxiter=10000, popsize=8)
print(r)
x = r.x
model = teqp.make_model(ppo.build_JSON(x))

Tc, rhoc = model.solve_pure_critical(400, 5000)
print(Tc, rhoc)
anc = teqp.build_ancillaries(model, Tc, rhoc, 0.5*Tc)

def _get_SOS(model, T, rho, z, *, R, M, Ao20):
    """ Helper function to calculate speed of sound """
    Ar0n = model.get_Ar02n(T, rho, z)
    Ar01 = Ar0n[1]; Ar02 = Ar0n[2]
    Ar11 = model.get_Ar11(T, rho, z)
    Ar20 = model.get_Ar20(T, rho, z)

    # M*w^2/(R*T) where w is the speed of sound
    # from the definition w = sqrt(dp/drho|s)
    Mw2RT = 1 + 2*Ar01 + Ar02 - (1 + Ar01 - Ar11)**2/(Ao20 + Ar20)
    if Mw2RT < 0:
        return 1e6
    w = (Mw2RT*R*T/M)**0.5
    return w

TcREF = CP.PropsSI('Tcrit', FLD)
Tsverify = np.linspace(100, TcREF*0.99999, 1000)
RHOL, RHOV, PPP, WWW = [],[],[],[]
z = np.array([1.0])
for T in Tsverify:
    rhoL, rhoV = model.pure_VLE_T(T, anc.rhoL(T)*1.01, anc.rhoV(T)*0.9, 10)
    pL = rhoL*model.get_R(z)*T*(1+model.get_Ar01(T, rhoL, z))
    pV = rhoV*model.get_R(z)*T*(1+model.get_Ar01(T, rhoV, z))
    RHOL.append(rhoL)
    RHOV.append(rhoV)
    PPP.append(pL)
    AS = CP.AbstractState('HEOS',FLD)
    AS.update(CP.DmolarT_INPUTS, 1e-10, T)
    Ao20 = AS.d2alpha0_dTau2()*AS.tau()**2
    WWW.append(_get_SOS(model, T, rhoL, z, R=8.314462618,M=CP.PropsSI('molemass',FLD),Ao20=Ao20))

# Plot the T-rho for VLE
line, = plt.plot(np.array(RHOL), Tsverify)
plt.plot(np.array(RHOV), Tsverify, color=line.get_color())
for Q in [0, 1]:
    D = CP.PropsSI('Dmolar','T',Tsverify,'Q',Q,FLD)
    plt.plot(D, Tsverify, lw=2, color='k')
plt.gca().set(xlabel=r'$\rho$ / mol/m$^3$', ylabel='$T$ / K')

# And a deviation plot, much closer to the critical point
# than the fitting region
fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(10,5), sharex=True)

ax1.plot(Tsverify, (np.array(PPP)/CP.PropsSI('P','T',Tsverify,'Q',0,FLD)-1)*100)
ax1.set(ylabel=r'$(p_{fit}/p_{\rm pexp}-1)\times 100$')

ax2.plot(Tsverify, (np.array(RHOL)/CP.PropsSI('Dmolar','T',Tsverify,'Q',0,FLD)-1)*100)
ax2.set(ylabel=r'$(\rho_{fit}/\rho_{\rm pexp}-1)\times 100$')

ax3.plot(Tsverify, (np.array(WWW)/CP.PropsSI('speed_of_sound','T',Tsverify,'Q',0,FLD)-1)*100)
ax3.set(ylabel=r'$(w_{fit}/w_{\rm pexp}-1)\times 100$', xlabel='$T$ / K')

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[2], line 2
      1 FLD = 'Methane'
----> 2 ppo = teqp.paramopt.PureParameterOptimizer(template, pointers)
      4 Ts = np.linspace(100, 170, 10)
      6 # Generate some artificial pseudo-experimental data from 
      7 # the reference EOS as implemented in CoolProp

AttributeError: module 'teqp' has no attribute 'paramopt'

Fit Pure Fluid Parameters¶

Related Topics

This Page