在 NumPy 中确定静力平衡

在分析物理结构时，了解保持其稳定的力学原理至关重要。作用在楼板、梁或任何其他结构上的力会产生反作用力和力矩。这些反作用力是结构抵抗运动而不破坏的。在结构即使受到力作用也不会移动的情况下，牛顿第二定律指出，系统中的加速度和所有方向上的合力都必须为零。您可以使用 NumPy 数组来表示和求解此概念。

您将做什么：

• 在本教程中，您将使用 NumPy 创建向量和力矩，使用 NumPy 数组。

• 解决涉及电缆和地板支撑结构的问题。

• 编写 NumPy 矩阵以隔离未知数。

• 使用 NumPy 函数执行线性代数运算。

您将学到什么：

• 如何使用 NumPy 表示点、向量和力矩。

• 如何找到向量的法向量

• 使用 NumPy 计算矩阵计算

您需要什么：

• NumPy

• Matplotlib

使用以下命令导入

import numpy as np
import matplotlib.pyplot as plt

在本教程中，您将使用以下 NumPy 工具

• np.linalg.norm：此函数确定向量大小的度量

• np.cross：此函数接受两个矩阵并产生叉积

使用牛顿第二定律求解平衡

您的模型由受力矩之和作用的梁组成。您可以使用牛顿第二定律开始分析此系统

\[\sum{\text{力}} = \text{质量} \times \text{加速度}.\]

为了简化所看的示例，假设它们是静态的，加速度 \(=0\)。由于我们的系统存在于三维空间中，因此应考虑在每个维度上施加的力。这意味着您可以将这些力表示为向量。对于力矩，您也会得出相同的结论，力矩是由作用在物体质心一定距离处的力产生的。

假设力 \(F\) 表示为三维向量

\[F = (F_x, F_y, F_z)\]

其中每个分量表示在每个相应方向上施加的力的幅度。还假设向量

\[r = (r_x, r_y, r_z)\]

中的每个分量都是力每个分量作用点与系统质心的距离。然后，力矩可以通过以下方式计算

\[r \times F = (r_x, r_y, r_z) \times (F_x, F_y, F_z).\]

从一些简单的力向量示例开始

forceA = np.array([1, 0, 0])
forceB = np.array([0, 1, 0])
print("Force A =", forceA)
print("Force B =", forceB)

Force A = [1 0 0]
Force B = [0 1 0]

这将 forceA 定义为幅度为 1 的 x 方向向量，并将 forceB 定义为幅度为 1 的 y 方向向量。

可视化这些力可能有助于更好地理解它们是如何相互作用的。Matplotlib 是一个具有可视化工具的库，可用于此目的。风筝图将用于演示三维向量，但该库也可用于二维演示。

fig = plt.figure()

d3 = fig.add_subplot(projection="3d")

d3.set_xlim(-1, 1)
d3.set_ylim(-1, 1)
d3.set_zlim(-1, 1)

x, y, z = np.array([0, 0, 0])  # defining the point of application.  Make it the origin

u, v, w = forceA  # breaking the force vector into individual components
d3.quiver(x, y, z, u, v, w, color="r", label="forceA")

u, v, w = forceB
d3.quiver(x, y, z, u, v, w, color="b", label="forceB")

plt.legend()
plt.show()

有两个力从一个点发出。为了简化此问题，您可以将它们加在一起以找到合力。请注意，forceA 和 forceB 都是三维向量，由具有三个分量的 NumPy 数组表示。由于 NumPy 旨在简化和优化向量之间的运算，因此您可以轻松计算这两个向量的和，如下所示

forceC = forceA + forceB
print("Force C =", forceC)

Force C = [1 1 0]

力 C 现在充当表示 A 和 B 的单个力。您可以绘制它以查看结果。

fig = plt.figure()

d3 = fig.add_subplot(projection="3d")

d3.set_xlim(-1, 1)
d3.set_ylim(-1, 1)
d3.set_zlim(-1, 1)

x, y, z = np.array([0, 0, 0])

u, v, w = forceA
d3.quiver(x, y, z, u, v, w, color="r", label="forceA")
u, v, w = forceB
d3.quiver(x, y, z, u, v, w, color="b", label="forceB")
u, v, w = forceC
d3.quiver(x, y, z, u, v, w, color="g", label="forceC")

plt.legend()
plt.show()

但是，目标是平衡。这意味着您希望合力为 \( (0, 0, 0) \)，否则您的物体将经历加速度。因此，需要另一个力来抵消先前的力。

您可以将此问题写为 \(A+B+R=0\)，其中 \(R\) 是解决问题的反作用力。

在此示例中，这意味着

\[(1, 0, 0) + (0, 1, 0) + (R_x, R_y, R_z) = (0, 0, 0)\]

分解为 x、y 和 z 分量，这将给您

\[\begin{split}\begin{cases} 1+0+R_x=0\\ 0+1+R_y=0\\ 0+0+R_z=0 \end{cases}\end{split}\]

求解 \(R_x\)、\(R_y\) 和 \(R_z\) 将给您一个向量 \(R\) 为 \( (-1, -1, 0) \)。

如果绘制，先前示例中看到的力应被抵消。只有当没有剩余的力时，系统才被认为处于平衡状态。

R = np.array([-1, -1, 0])

fig = plt.figure()

d3.set_xlim(-1, 1)
d3.set_ylim(-1, 1)
d3.set_zlim(-1, 1)

d3 = fig.add_subplot(projection="3d")

x, y, z = np.array([0, 0, 0])

u, v, w = forceA + forceB + R  # add them all together for sum of forces
d3.quiver(x, y, z, u, v, w)

plt.show()

空的图表表示没有外力。这表示系统处于平衡状态。

将平衡求解为力矩之和

接下来让我们看一个更复杂的应用。当力不都作用在同一点时，就会产生力矩。

类似于力，这些力矩也必须总和为零，否则将经历旋转加速度。类似于力的总和，这会为空间中的每个坐标方向创建线性方程。

一个简单的例子是作用在固定在地面上的静止杆上的力。杆不会移动，因此它必须施加反作用力。杆也不会旋转，因此它还必须产生反作用力矩。求解反作用力和力矩。

假设垂直于杆基部的 5N 力作用在杆上 2m 以上。

f = 5  # Force in newtons
L = 2  # Length of the pole

R = 0 - f
M = 0 - f * L
print("Reaction force =", R)
print("Reaction moment =", M)

Reaction force = -5
Reaction moment = -10

使用物理特性查找值

假设力不是垂直于梁作用，而是通过也连接到地面的电线作用到我们的杆上。给定该绳索的张力，您只需使用这些物体的物理位置即可解决此问题。

响应作用在杆上的力，基部产生了 x 和 y 方向的反作用力，以及一个反作用力矩。

将杆的底部表示为原点。现在，假设绳索连接到地面 3m 的 x 方向和杆上 2m 的 z 方向。

将这些点在空间中定义为 NumPy 数组，然后使用这些数组来查找方向向量。

poleBase = np.array([0, 0, 0])
cordBase = np.array([3, 0, 0])
cordConnection = np.array([0, 0, 2])

poleDirection = cordConnection - poleBase
print("Pole direction =", poleDirection)
cordDirection = cordBase - cordConnection
print("Cord direction =", cordDirection)

Pole direction = [0 0 2]
Cord direction = [ 3  0 -2]

为了将这些向量与力相关联，您需要将它们转换为单位向量。单位向量的大小为 1，并且只传达力的方向。

cordUnit = cordDirection / np.linalg.norm(cordDirection)
print("Cord unit vector =", cordUnit)

Cord unit vector = [ 0.83205029  0.         -0.5547002 ]

然后，您可以将此方向乘以力的幅度以找到力向量。

假设绳索的张力为 5N

cordTension = 5
forceCord = cordUnit * cordTension
print("Force from the cord =", forceCord)

Force from the cord = [ 4.16025147  0.         -2.77350098]

为了找到力矩，您需要力向量和半径的叉积。

momentCord = np.cross(forceCord, poleDirection)
print("Moment from the cord =", momentCord)

Moment from the cord = [ 0.         -8.32050294  0.        ]

现在您需要做的就是找到反作用力和力矩。

equilibrium = np.array([0, 0, 0])
R = equilibrium - forceCord
M = equilibrium - momentCord
print("Reaction force =", R)
print("Reaction moment =", M)

Reaction force = [-4.16025147  0.          2.77350098]
Reaction moment = [0.         8.32050294 0.        ]

另一个例子

让我们看一个稍微复杂一点的模型。在此示例中，您将观察带有两个电缆和一个施加力的梁。这次，您需要找到绳索的张力和梁的反作用力。(来源：工程师矢量力学：静力学，问题 4.106)

将距离 a 定义为 3 米

像以前一样，首先将每个相关点的坐标定义为一个数组。

A = np.array([0, 0, 0])
B = np.array([0, 3, 0])
C = np.array([0, 6, 0])
D = np.array([1.5, 0, -3])
E = np.array([1.5, 0, 3])
F = np.array([-3, 0, 2])

根据这些方程，您首先使用单位向量确定向量方向。

AB = B - C
AC = C - A
BD = D - B
BE = E - B
CF = F - C

UnitBD = BD / np.linalg.norm(BD)
UnitBE = BE / np.linalg.norm(BE)
UnitCF = CF / np.linalg.norm(CF)

RadBD = np.cross(AB, UnitBD)
RadBE = np.cross(AB, UnitBE)
RadCF = np.cross(AC, UnitCF)

这使您可以将张力 (T) 和反作用力 (R) 表示为

\[\begin{split}\left[ \begin{array} ~1/3 & 1/3 & 1 & 0 & 0\\ -2/3 & -2/3 & 0 & 1 & 0\\ -2/3 & 2/3 & 0 & 0 & 1\\ \end{array} \right] \left[ \begin{array} ~T_{BD}\\ T_{BE}\\ R_{x}\\ R_{y}\\ R_{z}\\ \end{array} \right] = \left[ \begin{array} ~195\\ 390\\ -130\\ \end{array} \right]\end{split}\]

以及力矩为

\[\begin{split}\left[ \begin{array} ~2 & -2\\ 1 & 1\\ \end{array} \right] \left[ \begin{array} ~T_{BD}\\ T_{BE}\\ \end{array} \right] = \left[ \begin{array} ~780\\ 1170\\ \end{array} \right]\end{split}\]

其中 \(T\) 是相应绳索的张力，\(R\) 是相应方向的反作用力。然后您只需六个方程

\(\sum F_{x} = 0 = T_{BE}/3+T_{BD}/3-195+R_{x}\)

\(\sum F_{y} = 0 = (-\frac{2}{3})T_{BE}-\frac{2}{3}T_{BD}-390+R_{y}\)

\(\sum F_{z} = 0 = (-\frac{2}{3})T_{BE}+\frac{2}{3}T_{BD}+130+R_{z}\)

\(\sum M_{x} = 0 = 780+2T_{BE}-2T_{BD}\)

\(\sum M_{z} = 0 = 1170-T_{BE}-T_{BD}\)

您现在有五个未知数和五个方程，可以求解

\(\ T_{BD} = 780N\)

\(\ T_{BE} = 390N\)

\(\ R_{x} = -195N\)

\(\ R_{y} = 1170N\)

\(\ R_{z} = 130N\)

总结

您已经学习了如何使用数组来表示三维空间中的点、力矩和力。数组中的每个条目都可以用来表示分解为方向分量的物理特性。然后，可以使用 NumPy 函数轻松地操作这些特性。

其他应用

此过程可以应用于动力学问题或任何数量的维度。本教程中的示例假设了静力平衡中的三维问题。这些方法可以轻松地用于更不同的问题。更多或更少的维度需要更大的或更小的数组来表示。在经历加速度的系统中，速度和加速度也可以表示为向量。

参考资料

1. 工程师矢量力学：静力学 (Beer & Johnston & Mazurek)

2. NumPy 参考