From f55190ca6a2e24b93c404045657c189fc157bf3a Mon Sep 17 00:00:00 2001
From: Marcel Frommelt <lrt86919@stud.uni-stuttgart.de>
Date: Mon, 3 May 2021 17:41:17 +0900
Subject: [PATCH] including several drag-models the user can choose from

---
 models/drag.py | 73 +++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 67 insertions(+), 6 deletions(-)

diff --git a/models/drag.py b/models/drag.py
index 38a4b2c..bbc7438 100644
--- a/models/drag.py
+++ b/models/drag.py
@@ -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