DPD shear viscosity plot

In this file the script that generate interactive plot as well as the plot presented.

Plot of viscosity for different values of deccorelation length

In [32]:
import numpy as np
from numpy import *
import glob
import copy
from bokeh.plotting import figure, output_notebook, show, output_file
output_notebook()
def sort2ar(a,b):
    a1=copy.copy(a)
    return [x for (a,x) in sorted(zip(a,b))]

inputdir="/home/azarnykh/Documents/PRL_article/article/postproc/visc_fit/velocity/sec3_alpha_0.0/"
etas=glob.glob(inputdir+"average_acf_*/ETAS.dat")
etasexp=glob.glob(inputdir+"average_acf_*/ETAS_EXP.dat")
etasint=glob.glob(inputdir+"average_acf_*/ETAS_INT.dat")

templist=[]

for ii,etai in enumerate(etas):
    m=1.0
    rc=1.0
    nden=4.0
    sigi=3.0
    aij=float(etai.split("/")[-3].split("_")[-1])
    temp_om=float(etai.split("/")[-2].split("_")[-1])
    temp = pow((pi*m*nden*rc**4*sigi**2*sqrt(m)/(45.0*temp_om)),2.0/3.0)
    templist.append(temp)

etas=sort2ar(templist,etas)
etasexp=sort2ar(templist,etasexp)
etasint=sort2ar(templist,etasint)
p = figure(width=500, height=500, tools="pan,wheel_zoom,box_zoom,reset,previewsave", x_axis_type="log",y_axis_type="log", title="Zero conservative force DPD effective shear viscosity")
#p = figure(width=500, height=500)
trace=[]
coll=["purple","red","brown","orange","gray","blue","black","green"]
for ii,etai in enumerate(etas):
    m=1.0
    rc=1.0
    nden=4.0
    sigi=3.0
    aij=float(etai.split("/")[-3].split("_")[-1])
    temp=float(etai.split("/")[-2].split("_")[-1])
    SCALE_zz=45.0*temp/pi/nden/sigi**2/m/rc**3
    OMEGA_0=pi*m*nden*rc**4*sigi**2*sqrt(m/temp)/(45.0*temp)
    l0=45.0*temp**1.5/(pi*m**1.5*nden*rc**3*sigi**2)

    et0=np.loadtxt(etai)
    et1=np.loadtxt(etasexp[ii])
    et2=np.loadtxt(etasint[ii])

    xscale=l0
    yscale=SCALE_zz/l0**2/4.0
    a2=1.0/35.0
    est=OMEGA_0**2*a2+0.5
    qrc=2.0*pi
    ql0=2.0*pi/l0

    if OMEGA_0>1:
        lname=str(int(round(OMEGA_0,1)))
    else:
        lname=str(round(OMEGA_0,2))
    x = array(et0[:,0])*xscale
    y = et0[:,1]*yscale
    xmax1=et0[0,0]*xscale
    p.line((ql0*xscale,ql0*xscale),(1e-3,100),line_color="gray",line_width = 1)
    p.line((1e-3, xmax1), (est,est), line_dash=[4, 4],color=coll[ii])
    p.line(x, y,line_color=coll[ii],line_width=2,legend=r'DPD, '+lname)
    p.circle(x, y,line_color=coll[ii],fill_color=None, size=7, alpha=0.5)
Loading BokehJS ...

Plot Langevin effective shear viscosity

In [33]:
import scipy.integrate
import scipy
m=1.
kT=0.25
be=0.1
mfp=(1./be)*sqrt(kT)

def fitfunc11(t,k):
    return exp(\
               -k**2 * kT/(m * be**2) *( exp(-be*t)+t*be-1)\
               -be   * t
    )
 
fitfunc1=fitfunc11
def yf(k):
    func = lambda t: fitfunc1(t,k)
    I=scipy.integrate.quad(func, 0., Inf)
    return I

kk=arange(0.001,0.1,0.001)
kk2=arange(0.1,1.,0.1)
kk3=arange(1.,100.,0.1)
kk=np.append(kk,kk2)
kk=np.append(kk,kk3)
ypl=[]
ypl2=[]
for i in kk:
    #print i
    a1=yf(i)[0]
    ypl.append(a1)
    l0=sqrt(kT)*1./be
kk=kk[kk<201]
Kn=kk
knvisc=(0.619/(0.619+(Kn))*(1.0+0.31*(Kn)/(0.785+1.152*(Kn)+(Kn)**2)))


l0=sqrt(kT)*1./be
t0=1./be
yle=t0/l0**2/np.array(ypl)/kk**2
xle=kk*l0
p.line(xle[yle<201], yle[yle<201],line_color="black",line_width=2,legend=r'LE theory')
Out[33]:
GlyphRenderer(
id = '5b3c3b33-9ffc-4404-aea3-0a64f249f2f7', …)

Plot scaling laws

In [34]:
Kn1=Kn[:120]
Kn2=Kn[120:]
yl1=5*Kn1**(-2)
yl2=1.25*Kn2**(-1)
p.line(Kn1[yl1<201], yl1[yl1<201], line_dash=[4, 4], color='gray')
p.line(Kn2[yl2<201], yl2[yl2<201], line_dash=[4, 4], color='gray')
Out[34]:
GlyphRenderer(
id = '82751484-2057-4bd7-888e-698aa7fbf80b', …)

Annotation

In [35]:
from bokeh.models import Label
text1=Label(x=2, y=0.01,text="q_l0")
p.add_layout(text1)
from bokeh.models import Label
text1=Label(x=1, y=20,text="q^-2")
p.add_layout(text1)
from bokeh.models import Label
text1=Label(x=100, y=0.02,text="q^-1")
p.add_layout(text1)

Output

In [36]:
show(p)