Merge pull request #314 from ls1mardyn/improvedPythonScripts

Revision of some python scripts in ppls1

tools/ppls1/ppls1/exp/chp.py

@@ -85,12 +85,8 @@ def exp_chp_bin_DF(fname, chp, append=False):
         writeMode = 'wb'
 
     with open(fname, writeMode) as f:
-        ba=bytearray()
-        for pi, row in chp.iterrows():
-            ba.extend(pack('<QIddddddddddddd',chp['pid'][pi],chp['cid'][pi],chp['rx'][pi],chp['ry'][pi],chp['rz'][pi],
-                                              chp['vx'][pi], chp['vy'][pi], chp['vz'][pi],chp['q0'][pi],chp['q1'][pi],
-                                              chp['q2'][pi], chp['q3'][pi], chp['Dx'][pi],chp['Dy'][pi],chp['Dz'][pi] ))
-        f.write(ba)
+        records = chp.to_records(index=False)
+        records.tofile(f)
         f.close()
 
 #%% Export ASCII ls1 checkpoint (representation: data frame)

tools/ppls1/ppls1/fluids/shear_visco.py

@@ -136,6 +136,7 @@ def eta(T,rho,fluid):
 
     na=6.02214076e23
     kb=1.380649e-23
+    u_mass=1.660539e-27
 
     if fluid == 'LJTS':
         fluid = 'Argon'
@@ -164,8 +165,8 @@ def eta(T,rho,fluid):
         T = T/tc*150.687
         rho = rho/sig**3*1e30/na*1e-3
         # refTime, refLambda: tref, lref
-        tref=sig*1e-10*np.sqrt(mass*1.660538e-27/(kb*eps))
-        eta_ref=1e3*kb/(sig*1e-10*tref) # [kg/m-s]  TODO
+        tref=sig*1e-10*np.sqrt(mass*u_mass/(kb*eps))
+        eta_ref=1e6*mass*u_mass/(sig*1e-10*tref) # Lemmon gives [uPa s] (see Table V in Paper)
     elif units == 'SI':
         pass
     else:
@@ -199,12 +200,59 @@ def eta_lauten(T,rho):
             eta_l = eta_l + coeff[i,0]*(T**(coeff[i,1]))*(rho**(coeff[i,2]))
     return eta_l
 
+
+#%% Get dynamic viscosity of LJfull fluid with correlation by Galliero, Ind. Eng. Chem. Res., vol. 44, 2005
+def eta_galliero(T,rho):
+    '''
+    Get dynamic viscosity of LJfull fluid (Galliero)
+
+    :param float T: Temperature
+    :param float rho: Density
+    :return: float eta_l: dynamic viscosity
+    '''
+
+    def omega_22(T):
+        a_22=1.16145
+        b_22=0.14874
+        c_22=0.52487
+        d_22=0.7732
+        e_22=2.16178
+        f_22=2.43787
+        omega_22 = a_22/T**b_22+c_22*np.exp(-d_22*T)+e_22*np.exp(-f_22*T)
+        return omega_22
+
+    b1=0.062692
+    b2=4.095577
+    b3=-8.743269e-6
+    b4=11.12492
+    b5=2.542477e-6
+    b6=14.863984
+
+    tmp1 = 0.17630924*np.sqrt(T)/omega_22(T)
+    tmp2 = b1*(np.exp(b2*rho)-1)+b3*(np.exp(b4*rho)-1)+b5*(np.exp(b6*rho)-1)/T**2
+
+    return tmp1+tmp2
+
+
+#%% Tests
 if __name__ == '__main__':
-    print('Running test with LJTS ...')
     fluid = 'LJTS'
-    T = 0.8
-    rho = 0.55
-
+    T = 0.95
+    rho = 0.7
+    print(f'Running test with {fluid} ...')
     print('Eta Lemmon:         '+str(eta_lemmon(T,rho,fluid)))
     print('Eta Lautenschlager: '+str(eta_lauten(T,rho)))
 
+    fluid = 'LJfull'
+    T = 1.1
+    rho = 0.7
+    print(f'Running test with {fluid} ...')
+    print('Eta Lemmon:   '+str(eta_lemmon(T,rho,fluid)))
+    print('Eta Galliero: '+str(eta_galliero(T,rho)))
+
+    fluid = 'Argon'
+    T = 200.0
+    rho = 10.0
+    print(f'Running test with {fluid} ...')
+    print('Eta Lemmon:     '+str(eta_lemmon(T,rho,fluid)))
+    print('Eta Lemmon Lit: 25.5662')  # Value from Table V in Lemmon2004

tools/ppls1/ppls1/fluids/therm_cond.py

@@ -139,6 +139,7 @@ def lamr(tau,delta,fluid):
 
     na=6.02214076e23
     kb=1.380649e-23
+    u_mass=1.660539e-27
 
     if fluid == 'LJTS':
         fluid = 'Argon'
@@ -167,8 +168,8 @@ def lamr(tau,delta,fluid):
         T = T/tc*150.687
         rho = rho/sig**3*1e30/na*1e-3
         # refTime, refLambda: tref, lref
-        tref=sig*1e-10*np.sqrt(mass*1.660538e-27/(kb*eps))
-        lref=1e3*kb/(sig*1e-10*tref) # [mW/m-K]
+        tref=sig*1e-10*np.sqrt(mass*u_mass/(kb*eps))
+        lref=1e3*kb/(sig*1e-10*tref) # Lemmon gives [mW/m-K] (see Table V in Paper)
     elif units == 'SI':
         pass
     else:
@@ -254,36 +255,19 @@ def lambda_fernandez(T,rho,L,QQ):
             lambda_f = lambda_f + coeff[i,j]*(T**(i))*(rho**(j))
     return lambda_f
 
-#%% Get thermal conductivity of LJTS fluid with data/fit by Guevara/Homes
-def lambda_guevara_homes(T,rho):
-    '''
-    Get thermal conductivity of LJTS fluid (Guevara/Homes)
-
-    :param float T: Temperature
-    :param float rho: Density
-    :return: float lambda_l: Thermal conductivity
-    '''
-
-    p00=-364.1
-    p10=-1003
-    p01=2292
-    p20=-21.11
-    p11=2238
-    p02=-3751
-    p30=601.2
-    p21=-1482
-    p12=13.63
-    p03=1494
-    lambda_gh = p00 + p10*T + p01*rho + p20*T**2 + p11*T*rho + p02*rho**2 + p30*T**3 + p21*T**2*rho + p12*T*rho**2 + p03*rho**3
-    return lambda_gh
 
+#%% Tests
 if __name__ == '__main__':
-    print('Running test with LJTS ...')
     fluid = 'LJTS'
-    T = 0.8
-    rho = 0.55
-
+    T = 0.9
+    rho = 0.75
+    print(f'Running test with {fluid} ...')
     print('Lambda Lemmon:         '+str(lambda_lemmon(T,rho,fluid)))
     print('Lambda Lautenschlager: '+str(lambda_lauten(T,rho)))
-    print('Lambda Guevara/Homes:  '+str(lambda_guevara_homes(T,rho)))
 
+    fluid = 'Argon'
+    T = 100.0
+    rho = 33.0
+    print(f'Running test with {fluid} ...')
+    print('Eta Lemmon:     '+str(lambda_lemmon(T,rho,fluid)))
+    print('Eta Lemmon Lit: 111.266')  # Value from Table V in Lemmon2004