最初，我有以下等式：
2mgv×sin(\alpha) = CdA×\rho(v^2 + v_{wind}^2 + 2vv_{wind}cos(\phi))^(3/2)
我可以将其表示为以下非线性方程：
(K × v)^(2/3) = v^2 + v_{wind} + 2vv_{wind}cos(\phi)
为了解决这个问题，我需要使用数值方法。
我尝试使用 scipy.optimize 中的 fsolve 在 Python 中为此编写代码。但是，我得到的结果不太理想。我还应该尝试什么？我应该使用不同的方法/包还是只是我的代码需要改进？我还发现结果高度依赖于 v_initial_guess。
请注意，我认为自己是编程初学者。
我也尝试编写代码，使用牛顿-拉夫逊法求解 v 方程，但不太成功。
这是我的代码：
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
编辑：这是我的代码现在的样子：
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)
9
最佳答案
4
您将无法使用牛顿-拉夫森公式（或其他任何公式）获得解决方案，因为…根本没有解决方案。
绘制原始方程式的 LHS 和 RHS 的图形：它们不相交（见下文）。因此，没有解。将 LHS 乘以 2，或者更自然地将 RHS 乘以 1/2，都不会改变这一点。
请提供原始方程的来源，因为从物理上讲，它根本不正确。你应该平衡力，而不是功率。如果你想平衡功率，那么将两个方程乘以实际速度 (v)，而不是左边的 v，右边的 v-v_wind。此外，如果 v 不是水平的，那么你的阻力无论如何都是错误的。
import numpy as np
import matplotlib.pyplot as plt
m = 80 # mass
g = 9.81 # gravity
alpha = np.radians(10) # incline
CdA = 0.65 # drag coefficient times area
rho = 1.225 # density of air
v_w = 10 # wind speed
phi = np.radians(30) # wind angle with the direction of motion
sin_alpha = np.sin(alpha)
cos_phi = np.cos(phi)
v = np.linspace( 0, 20, 101 )
LHS = m * g * v * sin_alpha # or twice this
#LHS = 2 * m * g * v * sin_alpha #
RHS = CdA * rho * ( v ** 2 + v_w ** 2 + 2 * v * v_w * cos_phi ) ** 1.5
plt.plot( v, LHS, label="LHS" )
plt.plot( v, RHS, label="RHS" )
plt.xlabel( "v" )
plt.ylabel( "power" )
plt.legend()
plt.show()
输出：
你没有告诉我们你的方程式从何而来。如果没有答案，我只能猜测了。
假设有一个物体，质量为 m，横截面积为 A，在水平横流 v w中以恒定速度 v 下降（水平面以下的 alpha 角） 。图中，红色表示力，蓝色表示速度。
如果速度稳定，则重量会平衡阻力。要实现此方向，则相对速度 U 必须沿所示方向（垂直向下）。
这给出了相对速度
如果知道风速，则速度和下降角度可由毕达哥拉斯定理和三角学给出：
另一方面，如果你知道下降角度，并且想找到风速和速度，那么
如果这一切都是错误的那么……请告诉我们你的方程式来自哪里！
|
您几乎已经弄清楚了情况。顺便说一句，我认为您在脚本中遗漏了 $\rho$。
但在转到物理部分之前，我插入了一个表示我的问题的图。在我看来，这个阻力曲线应该看起来有所不同。
因此，我思考物体在斜坡（\alpha$ 倾斜）上向上移动的方式在表达力量时如下：
$$P = P_{重力} + P_{摩擦力} + P_{空气阻力}$$
如果给定 $\alpha$，我想要找到 $v$ 速度，使得 P_{grav} = P_{air}。
$2mgvsin(\alpha) = CdA \rho × v_{rel}^3$
其中 $v_{rel} = \sqrt{v^2 + v_{wind}^2 + 2vv_{wind}×cos(\alpha)}$
替换后：
$(\frac {2mgvsin(\alpha)} {CdA\rho})^(2/3) = v^2 + v_{wind}^2 + 2vv_{wind}×cos(\phi)$
令 $K = \frac {2mgsin(\alpha)} {CdA \rho}
所以，
$$(Kv)^(2/3) = v^2 + v_{wind}^2 + 2vv_{wind}×cos(\phi)$$
正如我在第一篇帖子中看到的那样。
6
-
我不知道为什么我的 Tex 命令没有按应有的方式显示，抱歉。@lastchance
–
-
（我已在图形脚本中更正了 rho – 它不会改变结论）。您应该平衡力，而不是功率。功率 = 力 x 速度，而不是力 x（相对速度）。因此，CdA.rho.vrel^2.v，而不是 CdA.rho.vrel^3。（然后 v 取消 – 您最终平衡了力）
– -
@lastchance 我明白了，你说得对。我也喜欢你的方法，它似乎更简单。所以如果我理解正确的话，假设物体向上移动，U = \sqrt(( 2mg ×sin(alpha))/(CdA×rho)) 将是相对速度的解。而 v 的公式保持不变，对吗？
– -
我看到的是某个物体在空中下落，而不是在斜坡上移动。您的 vrel 公式似乎适用于水平面上的运动（并且无法平衡重量）。
– -
@lastchance 请查看上面的我的答案
–
|
$v_{rel} = \sqrt{v^2 + v_{wind}^2 + 2vv_{wind}×cos(\phi)}$ 是合成速度，定义了相对速度的大小。
$(v + v_{wind} * cos_phi)$ 这个因子代表行进方向的分量。
公式请参阅本文附录：
Sundström, D., Carlsson, P., Ståhl, F., & Tinnsten, M. (2013)。越野滑雪中步速策略的数值优化。结构与多学科优化，47，943-950。
|
谢谢你的帮助，@lastchance，它确实帮助我找到了我想要的解决方案。
请查看我的代码，该代码用于查找骑行中的翻转坡度（和速度），在此之后重力功率大于阻力功率（或反之亦然），因此，就性能而言，归一化为重量（W/kg）的功率比 W/CdA“更重要”（或反之亦然）。我为这种情况编写了另一个代码，当给定坡度并且需要翻转功率（或相应的速度）时。
在第一个工作版本中我仍然使用 fsolve，但后来我决定通过使用割线法来改变为一种“更好”的方法。
import numpy as np
import matplotlib.pyplot as plt
#what is the turnover gradient for a given PO (speed)?
P = 360
m = 80
g = 9.81
mu = 0.005
CdA = 0.321
rho = 1.225
v_wind = 4
phi = np.radians(60)
cos_phi = np.cos(phi)
def power_equation(v, v_wind):
alpha = np.arcsin((CdA * rho * (v + v_wind*cos_phi)**2) / (2 * m * g))
gravity = m * g * np.sin(alpha) * v
friction = mu * m * g * np.cos(alpha) * v
air = 0.5 * CdA * rho * (v + v_wind * np.cos(phi))**2 * v
return gravity + friction + air - P
def secant_method(func, x0, x1, args=(), tol=1e-8, max_iter=1000):
for i in range(max_iter):
f0 = func(x0, *args)
f1 = func(x1, *args)
if abs(f1 - f0) < tol:
return x1
x_next = x1 - f1 * (x1 - x0) / (f1 - f0)
x0 = x1
x1 = x_next
raise RuntimeError("Secant method did not converge")
x0 = 10.0
x1 = 11.0
v_solution = secant_method(power_equation, x0, x1, args=(v_wind,))
print(f"\nP [W]= {P} \nv_turnover [m/s] = {v_solution:.2f}\nabove which drag power (W/CdA) matters more")
alpha1 = np.arcsin((CdA * rho * (v_solution + v_wind*cos_phi)**2) / (2 * m * g))
alpha1_degrees = np.degrees(alpha1)
alpha1_percentage = np.tan(alpha1) * 100
print(f"turnover gradient [deg]: {alpha1_degrees}")
print(f"turnover gradient [%]: {alpha1_percentage}\n above which gravity power (W/kg) matters more")
P_grav = m * g * np.sin(alpha1) * v_solution
P_fric = mu * m * g * np.cos(alpha1) * v_solution
P_air = 0.5 * CdA * rho * (v_solution + v_wind * np.cos(phi))**2 * v_solution
print(f"\nP_gravity: {P_grav:.2f}")
print(f"P_friction: {P_fric:.2f}")
print(f"P_air: {P_air:.2f}")
total_power = P_grav + P_fric + P_air
print(f"Total Power: {total_power} W")
print(f"{(P_grav/total_power)*100:.2f} %, {(P_air/total_power)*100:.2f} %, {(P_fric/total_power)*100:.2f} %")
#PLOT
alpha_degrees = np.linspace(0, 45, 100)
alpha_values = np.radians(alpha_degrees)
velocity_values = []
gravity_power_values = []
drag_power_values = []
friction_power_values = []
for alpha in alpha_values:
def power_eq_with_alpha(v):
gravity = m * g * np.sin(alpha) * v
friction = mu * m * g * np.cos(alpha) * v
air = 0.5 * CdA * rho * (v + v_wind * cos_phi)**2 * v
return gravity + friction + air - P
x0 = 10.0
x1 = 11.0
try:
v_solution_array = secant_method(power_eq_with_alpha, x0, x1)
velocity_values.append(v_solution_array)
gravity_power_output = m * g * np.sin(alpha) * v_solution_array
friction_power_output = mu * m * g * np.cos(alpha) * v_solution_array
drag_power_output = 0.5 * CdA * rho * (v_solution_array + v_wind * cos_phi)**2 * v_solution_array
gravity_power_values.append(gravity_power_output)
friction_power_values.append(friction_power_output)
drag_power_values.append(drag_power_output)
except RuntimeError:
velocity_values.append(np.nan)
gravity_power_values.append(np.nan)
friction_power_values.append(np.nan)
drag_power_values.append(np.nan)
velocity_values = np.array(velocity_values)
gravity_power_values = np.array(gravity_power_values)
friction_power_values = np.array(friction_power_values)
drag_power_values = np.array(drag_power_values)
total_power_values = gravity_power_values + friction_power_values + drag_power_values
fig, ax1 = plt.subplots(figsize=(15, 6))
ax1.plot(alpha_degrees, total_power_values, label='Power output')
ax1.plot(alpha_degrees, gravity_power_values, label='P_grav')
ax1.plot(alpha_degrees, friction_power_values, label='P_fric')
ax1.plot(alpha_degrees, drag_power_values, label='P_air')
ax1.axvline(x=alpha1_degrees, color='grey', linestyle='--')
ax1.set_xlabel('Gradient (degrees)')
ax1.set_ylabel('Power (W)')
ax1.legend(loc='upper left')
ax1.set_title('Gradient vs Power Components and Speed')
ax2 = ax1.twinx()
ax2.plot(alpha_degrees, velocity_values, linestyle='--', label='Speed', color='orange')
ax2.axhline(y=v_solution, color='blue', linestyle='-.')
ax2.set_ylabel('Speed (m/s)')
ax2.legend(loc='upper right')
plt.show()
1
-
@lastchance 再次感谢您的时间和耐心！我知道我越野了几次……
–
|
–
I wasnt too successful.结果究竟如何？–
–
–
–
|