import numpy as np
import holoviews as hv
hv.extension('bokeh', logo=False)
rhosw = 1.65 # 水の水中比重
g = 9.8 # 重力加速度
def rubeyw0(dd, nu=10**(-6)):
tmp = 36.0*nu**2/rhosw/g/dd**3
F = np.sqrt(2.0/3.0 + tmp) - np.sqrt(tmp) if dd < 0.001 else np.sqrt(2.0/3.0)
return np.sqrt(rhosw*g*dd)*F
def erfc(xi):
def erf(xi):
x = np.abs(xi)
t = 1/(1+0.5*x)
tt = t*np.exp( -x**2
-1.26551223*t**0
+1.00002368*t**1
+0.37409196*t**2
+0.09678418*t**3
-0.18628806*t**4
+0.27886807*t**5
-1.13520398*t**6
+1.48851587*t**7
-0.82215223*t**8
+0.17087277*t**9
)
return 1.0 - tt if xi >= 0.0 else tt-1.0
return 1.0 - erf(xi)
Kは河床形態による係数. 0.0018を使用
def ca(dd, taus, w0, K = 0.0018):
alfa = 0.14
Bs = 0.143
eta = 0.5
def omega(taus):
aprime = Bs/taus - 1.0/eta
c1 = 0.5 / np.sqrt(np.pi) * np.exp(-aprime**2)
c2 = 0.5 * erfc(aprime)
if c2 <= 0.0 :
rr = 0.0
else:
rr =taus/Bs*c1/c2+ taus/Bs/eta - 1.0
return rr
us = np.sqrt(rhosw*g*dd*taus)
r = K*( alfa/(rhosw+1)*us/w0*omega(taus)/taus - 1.0 )
return 0.0 if r < 0.0 else r
r = np.arange(-1.3,2.1,0.1)
taus = 10**r
d = 0.173/1000
w0 = rubeyw0(d)
us = np.sqrt(rhosw*g*d*taus)
aaa = [ ca(d, tt, w0) for tt in taus]
hv.Curve((w0/us, aaa)).options(logx=True, logy=True)
WARNING:param.CurvePlot01361: Logarithmic axis range encountered value less than or equal to zero, please supply explicit lower-bound to override default of 0.000.
h = 5.0
ib = 1.0/1000.0
def cfig(h,ib,d):
us = np.sqrt(g*h*ib)
w0 = rubeyw0(d)
taus = us**2/rhosw/g/d
kappa=0.4
zeta = w0/kappa/us
y = np.arange(0.05,1.0,0.01) * h
a = 0.05*h
x = (y-a)/(h-a)
CbyCa = (a*(h-y)/y/(h-a))**zeta
Ca = ca(d, taus, w0)
return CbyCa*Ca*10**6, x
g1 = hv.Curve(cfig(h,ib,0.1/1000),label='0.1mm',kdims=['C']).options(logx=True, logy=False)
g2 = hv.Curve(cfig(h,ib,0.2/1000),label='0.2mm',kdims=['C']).options(logx=True, logy=False)
g3 = hv.Curve(cfig(h,ib,0.5/1000),label='0.5mm',kdims=['C']).options(logx=True, logy=False)
g4 = hv.Curve(cfig(h,ib,1.0/1000),label='1.0mm',kdims=['C']).options(logx=True, logy=False)
g5 = hv.Curve(cfig(h,ib,2.0/1000),label='2.0mm',kdims=['C']).options(logx=True, logy=False)
f1 = (g1*g2*g3*g4*g5).options(xlabel='土砂濃度[ppm]', ylabel='相対水深',width=500, legend_position='bottom_left').redim.range(y=(0,1.0))
f1
hv.save(f1,'suspendProf.html')