Source code for fastwater.general.physics
# -*- coding: utf-8 -*-
"""
Functions to calculate various physical parameters used in fluid modelling.
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import numpy as np
# Non-Standard Python Dependencies
# Local Module Dependencies
from fastwater.general.geometry import unitVect, circSign2D
from fastwater.general.meshTools import getTriangleArea
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
grav = 9.81 # Mean gravitional acceleration at Earths surface m / s^2
wgs84_a = 6378137.0 # WGS84 elipsoid major axis
wgs84_b = 6356752.3 # WGS84 elipsoid minor axis
wgs84_gm = 9.7803253359 # WGS84 equatorial gravitional acceleration
wgs84_gn = 0.001931850400 # WGS84 gravity numerator constant
wgs84_gd = 0.006694384442 # WGS84 gravity denominator constant
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
#def g_wgs84(lat: float):
[docs]def g_wgs84(lat):
"""
Calculate gravity based on latitude on the WGS84 spheroid.
Parameters
----------
lat : float
Latitude in decimal degrees.
Returns
-------
g : float
Gravitational acceleration at latitude in m/s^2.
"""
g = wgs84_gm * ((1.0 + wgs84_gn * np.square(np.sin(np.radians(lat))))/
np.sqrt(1.0 - wgs84_gd * np.square(np.sin(np.radians(lat)))))
return g
[docs]def circulation2D(x, y, vel_x, vel_y):
"""
Calculate the circulation around an element define by the nodes at
locations ( x , y ).
Parameters
----------
x : list[ float ]
x coordinates of element nodes [ m ].
y : list[ float ]
y coordinate of element nodes [ m ].
vel_x : list[ float ]
x component of velocity at element nodes [ m/s ].
vel_y : list[ float ]
y component of velocity at element nodes [ m/s ].
Returns
-------
circ : float
fluid circulation around element.
"""
npts = len(x)
x.append(x[0])
y.append(y[0])
pts = np.array([x,y]).transpose()
vel_x.append(vel_x[0])
vel_y.append(vel_y[0])
vels = np.array([vel_x,vel_y]).transpose()
edgeLen = np.zeros(npts,)
edgeTan = np.zeros((npts,2))
edgeVel = np.zeros((npts,2))
circ = 0.0
if npts <3:
circ = None
else:
x = np.zeros(npts+1,)
y = np.zeros(npts+1,)
for ipt in range(npts):
edgeLen[ipt] = np.linalg.norm(pts[ipt+1]-pts[ipt])
edgeTan[ipt,:] = unitVect(pts[ipt],pts[ipt+1])
edgeVel[ipt,:] = (vels[ipt,:] + vels[ipt+1,:]) / 2.0
area = getTriangleArea(pts[0],pts[1],pts[2])
circ = np.sum(np.sum(edgeTan*edgeVel,1)*edgeLen)/area
# ensure anticlockwise integration around edges
csgn = circSign2D(edgeTan[0],edgeTan[1])
circ = circ * csgn
return circ
[docs]def kinetic_energy_density_2D(mesh,depth,velx,vely,rho):
# K.E. density = 1/2 m V^2 / A
# = 1/2 (rho * Vol) V^2 / A
# = 1/2 (rho * A * d) V^2 / A
# = 1/2 rho d V^2
#
d = np.sum(depth[mesh['elems']]*mesh['cWghts'],1);
u = np.sum(velx[mesh['elems']]*mesh['cWghts'],1);
v = np.sum(vely[mesh['elems']]*mesh['cWghts'],1);
v2 = np.square(u)+np.square(v);
ke_dens = 0.5*rho*d*v2;
return ke_dens
[docs]def potential_energy_density_2D(mesh,elev,rho):
# P.E. density = m g h / A
# = (rho * Vol) g h / A
# = (rho * A * h) g h / A
# = rho g h^2
#
h = np.sum(elev[mesh['elems']]*mesh['cWghts'],1);
pe_dens = rho*grav*np.square(h);
return pe_dens
[docs]def frictionCoeff(C100):
return (2.0/3.0)*C100
[docs]def chezy(C100):
"""
Calculate the Chezy friction coefficient from the C_100 drag coeficient.
Parameters
----------
C100 : float
C_100 drag coefficient (based on the velocity 100cm above bed).
Returns
-------
Cz : float
Chezy friction coefficient.
"""
Cf = frictionCoeff(C100)
Cz = np.sqrt(2.0*grav/Cf)
return Cz
[docs]def manning(C100,h):
"""
Calculate the Manning friction coefficient based on the C_100 drag
coefficient and the water depth.
Parameters
----------
C100 : float
Drag coeffieient at 100cm above the beed.
h : float
water depth [ m ].
Returns
-------
m : float
Manning friction coefficient.
"""
Cf = frictionCoeff(C100)
m = np.sqrt(Cf*np.power(np.abs(h),1.0/3.0)/(2.0 * grav))
return m
[docs]def strickler(C100,h):
"""
Calculate the Strickler friction coefficient based on the C_100 drag
coefficient and the water depth.
Parameters
----------
C100 : float
Drag coeffieient at 100cm above the beed.
h : float
water depth [ m ].
Returns
-------
S : float
Strickler friction coefficient.
"""
Cf = frictionCoeff(C100)
S = np.sqrt(2.0*grav/(Cf*np.power(np.abs(h),1.0/3.0)))
return S
[docs]def nikuradse(C100,h):
"""
Calculate the Nikuradse friction coefficient based on the C_100
drag coeffieient and the water depth.
Parameters
----------
C100 : float
Drag coeffieient at 100cm above the beed.
h : float
water depth [ m ].
Returns
-------
ks : float
Nikuradse friction coefficient.
"""
Cz = chezy(C100)
ks = 12.0*np.abs(h)*np.exp(-Cz/7.83)
return ks