Originally, I had the following equation:
2mgv×sin(\alpha) = CdA×\rho(v^2 + v_{wind}^2 + 2vv_{wind}cos(\phi))^(3/2)
which I could express as the following non-linear equation:
(K × v)^(2/3) = v^2 + v_{wind} + 2vv_{wind}cos(\phi)
To solve this, I need to use a numerical approach.
I tried writing a code for this in Python using fsolve from scipy.optimize. However, the results I got are not too promising. What else should I try? Should I use a different approach/package or just my code needs to be improved? I also experienced that the result is highly dependent on v_initial_guess.
Please note, that I would consider myself as a beginner in programming.
I also tried writing a code to solve the equation for v using the Newton-Raphson method but I wasnt too successful.
Here is my code:
import numpy as np
from scipy.optimize import fsolve
m = 80
g = 9.81
alpha = np.radians(10) #incline
CdA = 0.65
rho = 1.225
v_w = 10
phi = np.radians(30) #wind angle with the direction of motion
sin_alpha = np.sin(alpha)
cos_phi = np.cos(phi)
def equation(v):
K = ((m * g * sin_alpha) / ((CdA * rho))**(2/3))
return K * v**(2/3) - v**2 - 2*v*v_w*cos_phi - v_w**2
v_initial_guess = 30
v_solution = fsolve(equation, v_initial_guess, xtol=1e-3)
print("v:", v_solution[0])type here
EDIT: This is what my code looks like now:
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
m = 80
g = 9.81
alpha = np.radians(2) # incline
CdA = 0.321
rho = 1.22
v_w = 5
phi = np.radians(180) # wind angle with the direction of motion
sin_alpha = np.sin(alpha)
cos_phi = np.cos(phi)
def lhs(v):
return m * g * v * sin_alpha
def rhs(v):
return 0.5 * CdA * rho * (v**2 + v_w**2 + 2*v*v_w*cos_phi)**(3)
def difference(v):
return lhs(v) - rhs(v)
# fsolve to find the intersection
v_initial_guess = 8
v_intersection = fsolve(difference, v_initial_guess)[0]
v_values = np.linspace(0.1, 50, 500)
lhs_values = lhs(v_values)
rhs_values = rhs(v_values)
plt.figure(figsize=(10, 6))
plt.plot(v_values, lhs_values, label='$2mgv\\sin(\\alpha)$', color='blue')
plt.plot(v_values, rhs_values, label='$CdA\\rho(v^2 + v_{wind}^2 + 2vv_{wind}\\cos(\\phi))^{3/2}$', color='red')
plt.xlabel('Velocity (v)')
plt.xlim(0, 20)
plt.title('LHS and RHS vs. Velocity')
plt.legend()
plt.grid(True)
plt.ylim(0, 2000)
plt.show()
print(f"The intersection occurs at v = {v_intersection:.2f} m/s")
P_grav_check = m *g * sin_alpha * v_intersection
P_air_check = 0.5 * CdA * rho * (v_intersection ** 2 + v_w ** 2 + 2 * v_intersection * v_w * cos_phi) ** (3)
print(P_grav_check)
print(P_air_check)
The solution, using the secant method is the following:
import numpy as np
import matplotlib.pyplot as plt
m = 60
g = 9.81
alpha = np.radians(2) # incline
CdA = 0.321
rho = 1.225
v_wind = 4
phi = np.radians(60) # wind angle with the direction of motion
mu = 0.005
cos_phi = np.cos(phi)
sin_alpha = np.sin(alpha)
lhs = m * g * sin_alpha # gravitational force
def rhs(v):
return 0.5 * CdA * rho * (v + v_wind *cos_phi)**2 # drag force
def difference(v):
return lhs - rhs(v)
def secant_method(func, v0, v1, tolerance=1e-6, max_iterations=1000):
for _ in range(max_iterations):
f_v0 = func(v0)
f_v1 = func(v1)
v_new = v1 - f_v1 * (v1 - v0) / (f_v1 - f_v0)
if abs(v_new - v1) < tolerance:
return v_new
v0, v1 = v1, v_new
raise ValueError("Secant method did not converge")
v_initial_guess_1 = 15
v_initial_guess_2 = 25
v_intersection = secant_method(difference, v_initial_guess_1, v_initial_guess_2)
v_values = np.linspace(0.1, 50, 500)
rhs_values = rhs(v_values)
P_grav_check = m * g * sin_alpha * v_intersection
P_air_check = 0.5 * CdA * rho * (v_intersection + v_wind * cos_phi)**2 * v_intersection
P_friction = mu * m * g * np.cos(alpha) * v_intersection
alpha_degree = np.degrees(alpha)
alpha_percentage = np.tan(alpha) *100
print(f"\nopposing component of v_wind: ",v_wind *cos_phi, "m/s")
print(f"The intersection occurs at v = {v_intersection:.4f} m/s")
print(f" at a {alpha_degree} degrees incline which equals to {alpha_percentage:.2f} %")
print(f"\nP_grav: {P_grav_check:.2f} W")
print(f"P_air_check: {P_air_check:.2f} W")
print(f"P_friction: {P_friction:.2f} W")
print(f"P = {P_grav_check + P_air_check + P_friction}")
plt.figure(figsize=(10, 6))
plt.plot(v_values, [lhs]*len(v_values), label='$mg\\sin(\\alpha)$', color='blue')
plt.plot(v_values, rhs_values, label='$0.5CdA\\rho ({v + v_{wind}\\cos(\\phi)})^2$', color='red')
plt.xlabel('Velocity (v)')
plt.xlim(0, 20)
plt.title('Gravitational vs. Drag Force')
plt.legend()
plt.grid(True)
plt.ylim(0, 200)
plt.show()