esbo_ds/BASTET
esbo_ds/BASTET
10
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

including several drag-models the user can choose from

Browse Source
This commit is contained in:
Marcel Christian Frommelt 2021-05-03 17:41:17 +09:00
parent 2317c2907c
commit f55190ca6a
1 changed files with 67 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
models/drag.py
View File

@ -2,22 +2,83 @@ import numpy as np
from input.user_input import *
def Cd_model(fr, re, a_top):
k_cd = 19.8
k1 = 1.50 * 10 ** 7
k2 = 0.86 * 10 ** (-7)
a_top0 = np.pi / 4 * 1.383 * V_design ** (2/3) # a_top0 for zero-pressure shape
def cd_Palumbo(fr, re, a_top):
k_cd = 19.8 # +/- 1.2
k1 = 1.50 * 10 ** 7 # +/- 0.78
k2 = 0.86 * 10 ** (-7) # +/- 0.43
a_top0 = np.pi / 4 * 1.383 ** 2 * V_design ** (2/3) # a_top0 for zero-pressure shape
value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2))
if value > 1.6:
value = 1.6
elif value <= 0.1:
value = 0.1
return value
c_d = 0.8 # 0.47
def cd_PalumboLow(fr, re, a_top):
k_cd = 19.8 - 1.2
k1 = (1.50 - 0.78) * 10 ** 7
k2 = (0.86 - 0.43) * 10 ** (-7)
a_top0 = np.pi / 4 * 1.383 ** 2 * V_design ** (2/3) # a_top0 for zero-pressure shape
value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2))
if value > 1.6:
value = 1.6
elif value <= 0.1:
value = 0.1
return value
def cd_PalumboHigh(fr, re, a_top):
k_cd = 19.8 + 1.2
k1 = (1.50 + 0.78) * 10 ** 7
k2 = (0.86 + 0.43) * 10 ** (-7)
a_top0 = np.pi / 4 * 1.383 ** 2 * V_design ** (2/3) # a_top0 for zero-pressure shape
value = 0.2 * k_cd/fr * (k1/re + k2 * re) * a_top / (a_top0 ** (3/2))
if value > 1.6:
value = 1.6
elif value <= 0.1:
value = 0.4
return value
def drag(c_d, rho_air, A_proj, v):
return 0.5 * c_d * rho_air * A_proj * v ** 2
def cd_sphere(Re):
"This function computes the drag coefficient of a sphere as a function of the Reynolds number Re."
# Curve fitted after fig . A -56 in Evett and Liu: "Fluid Mechanics and Hydraulics"
from numpy import log10, array, polyval
if Re <= 0.0:
CD = 0.0
elif Re > 8.0e6:
CD = 0.2
elif Re > 0.0 and Re <= 0.5:
CD = 24.0 / Re
elif Re > 0.5 and Re <= 100.0:
p = array([4.22, -14.05, 34.87, 0.658])
CD = polyval(p, 1.0 / Re)
elif Re > 100.0 and Re <= 1.0e4:
p = array([-30.41, 43.72, -17.08, 2.41])
CD = polyval(p, 1.0 / log10(Re))
elif Re > 1.0e4 and Re <= 3.35e5:
p = array([-0.1584, 2.031, -8.472, 11.932])
CD = polyval(p, log10(Re))
elif Re > 3.35e5 and Re <= 5.0e5:
x1 = log10(Re / 4.5e5)
CD = 91.08 * x1 ** 4 + 0.0764
else:
p = array([-0.06338, 1.1905, -7.332, 14.93])
CD = polyval(p, log10(Re))
return CD