Commit d4d72fde authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

water viscosity calculation simplifications

parent bd44ba52

pysen/physics/water_model.py

+13 −18
Original line number Diff line number Diff line
@@ -7,7 +7,8 @@ H2O/D2O Model Parameters 
"""



from scipy.special import spherical_jn

from numpy import asarray, arange, ones_like, newaxis, sqrt, exp, log, interp

from numpy import ( asarray, arange, ones_like, newaxis, sqrt, exp, log,

                    interp, polyval )

from numpy import sum as npsum



# Model parameters, see Table
@@ -404,28 +405,22 @@ 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

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

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

    """

    A = 802.25336e-3

    a = 3.4741e-3

    b =-1.7413e-5

    c = 2.7719e-8

    gamma = 1.53026

    T0 = 225.334 # K

    dt = T-T0

    return A*(dt+a*dt**2+b*dt**3+c*dt**4)**(-gamma)



    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

    """

    A = 885.60402e-3

    a = 2.7990e-3

    b =-1.6342e-5

    c = 2.9067e-8

    gamma = 1.55255

    T0 = 231.832 # K

    dt = T-T0

    return A*(dt+a*dt**2+b*dt**3+c*dt**4)**(-gamma)
 
    return _water_viscosity(T-231.832, 885.60402e-3, 1.55255,

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

pysen/revision.py

+2 −2
Original line number Diff line number Diff line
@@ -3,8 +3,8 @@ PySEN revision module 
"""

import sys

__version__  = "2.0"

__release__  = ".2b2"

__date__     = "July 2, 2025"

__release__  = ".2b3"

__date__     = "July 3, 2025"



def version(full=False):

    "get pysen version number"

test/test_physics_water_model.py

+7 −6
Original line number Diff line number Diff line
@@ -60,12 +60,13 @@ class WaterModelTestCase(unittest.TestCase): 
    def test_viscosity(self):

        "test viscosity"

        C = 273.15

        self.assertAlmostEqual(wm.eta_h2o( 5.0+C)*1e3, 1.519, 3)

        self.assertAlmostEqual(wm.eta_h2o(20.0+C)*1e3, 1.002, 3)

        self.assertAlmostEqual(wm.eta_h2o(70.0+C)*1e3, 0.404, 3)

        self.assertAlmostEqual(wm.eta_d2o( 5.0+C)*1e3, 1.988, 3)

        self.assertAlmostEqual(wm.eta_d2o(20.0+C)*1e3, 1.251, 3)

        self.assertAlmostEqual(wm.eta_d2o(70.0+C)*1e3, 0.474, 3)

        mili=1e-3

        self.assertAlmostEqual(wm.eta_h2o( 5.0+C)/mili, 1.51893, 5)

        self.assertAlmostEqual(wm.eta_h2o(20.0+C)/mili, 1.00203, 5)

        self.assertAlmostEqual(wm.eta_h2o(70.0+C)/mili, 0.40415, 5)

        self.assertAlmostEqual(wm.eta_d2o( 5.0+C)/mili, 1.98777, 5)

        self.assertAlmostEqual(wm.eta_d2o(20.0+C)/mili, 1.25136, 5)

        self.assertAlmostEqual(wm.eta_d2o(70.0+C)/mili, 0.47435, 5)



    def test_teixeira(self):

        "test Teixeira model"