Commit 91c162f4 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

water model cleaning

parent d1a7b836

pysen/physics/water_model.py

+20 −21
Original line number Diff line number Diff line
@@ -404,28 +404,11 @@ def tx_model(t, Q, T, components='all'): 




# ===============================================================

# H2O and D2O viscosity parameterization

# by Cho et al.  J. Phys. Chem. B 103 (1999), 1991-1994

# H2O and D2O density and viscosity parameterization

#

# Density: 2nd degree polynomial fit to NIST data (5-40C)

# Viscosity Cho et al.  J. Phys. Chem. B 103 (1999), 1991-1994

# ===============================================================

def _water_viscosity(dt,A, gamma, *acoef):

    """

    returns A*(acoef[0]*dt**4+acoef[1]*dt**3+acoef[2]*dt**2+dt)**(-gamma)

    """

    xt = (polyval(*acoef,dt)*dt+1)*dt

    return A*xt**(-gamma)



def eta_h2o(T):

    """returns water viscosity in  Pa*s

    T in Kelvin

    """

    return _water_viscosity(T-225.334, 802.25336e-3, 1.53026,

                            (2.7719e-8, -1.7413e-5, 3.4741e-3))

def eta_d2o(T):

    """returns heavy water viscosity in  Pa*s

    T in Kelvin

    """

    return _water_viscosity(T-231.832, 885.60402e-3, 1.55255,

                            (2.9067e-8, -1.6342e-5, 2.7990e-3))



# simple fit

def rho_h2o(T):
@@ -439,3 +422,19 @@ def rho_d2o(T): 
    temp = T-T0C  # convert to Celsius

    assert 5.0<=temp<=40.0, "temperature must be between 5C and 40C"

    return 1.10517 + 1.42619e-4*temp - 6.80952e-6*temp**2



# ===============================================================

# note: 1 cP = 1 mPa*s = 10^-3 Pa*s

def _water_viscosity(dt,A, gamma, *acoef):

    "returns A*(acoef[0]*dt**4+acoef[1]*dt**3+acoef[2]*dt**2+dt)**(-gamma)"

    xt = (polyval(*acoef,dt)*dt+1)*dt

    return A*xt**(-gamma)



def eta_h2o(T):

    "returns water viscosity in  Pa*s,  T in Kelvin "

    return _water_viscosity(T-225.334, 802.25336e-3, 1.53026,

                            (2.7719e-8, -1.7413e-5, 3.4741e-3))

def eta_d2o(T):

    "returns heavy water viscosity in  Pa*s, T in Kelvin "

    return _water_viscosity(T-231.832, 885.60402e-3, 1.55255,

                            (2.9067e-8, -1.6342e-5, 2.7990e-3))

test/test_physics_water_model.py

+4 −4
Original line number Diff line number Diff line
@@ -75,10 +75,10 @@ class WaterModelTestCase(unittest.TestCase): 
                        25.0: 0.9970, 30.0: 0.9957, 35.0: 0.9940 }

        exp_rho_d2o = { 10.0: 1.1059, 15.0: 1.1058, 20.0: 1.1053,

                        25.0: 1.1045, 30.0: 1.1033} 

        for t in exp_rho_h2o:

            self.assertAlmostEqual(wm.rho_h2o(t+C), exp_rho_h2o[t], 4, f"rho_h2o at {t}C failed")

        for t in exp_rho_d2o:

            self.assertAlmostEqual(wm.rho_d2o(t+C), exp_rho_d2o[t], 4, f"rho_d2o at {t}C failed")

        for t,v in exp_rho_h2o.items():

            self.assertAlmostEqual(wm.rho_h2o(t+C), v, 4, f"rho_h2o at {t}C failed")

        for t,v in exp_rho_d2o.items():

            self.assertAlmostEqual(wm.rho_d2o(t+C), v, 4, f"rho_d2o at {t}C failed")



    def test_teixeira(self):

        "test Teixeira model"