液体力#

液体力学的精确模拟超出了 MuJoCo 的范围，并且对于我们旨在支持的应用来说速度太慢。尽管如此，我们提供了两种唯象模型，足以模拟飞行和游泳等行为。这些模型是无状态的，即没有为周围流体分配额外的状态，但它们能够捕捉刚体在流体介质中运动的显著特征。

通过将密度 (density) 和粘度 (viscosity) 属性设置为正值来启用这两种模型。这些参数对应于介质的密度 \(\rho\) 和粘度 \(\beta\)。

基于惯性的模型仅使用粘度和密度，从等效惯性箱体推断几何形状。
基于椭球体的模型更为精细，使用几何体的椭球体近似。除了介质的全局粘度和密度外，此模型还为每个相互作用的几何体提供了 5 个可调参数。

提示

如数值积分部分所述，隐式积分在存在与速度相关的力时显著提高了仿真稳定性。下面描述的两种液体力模型都具有此特性，因此在使用液体力时建议使用 implicit 或 implicitfast 积分器。这两种模型所需的解析导数已完全实现。

惯性模型#

在此模型中，为了液体力学目的，每个物体的形状被假定为等效惯性箱体，该箱体也可以被可视化。对于质量为 \(\mathcal{M}\) 和惯性矩阵为 \(\mathcal{I}\) 的物体，等效惯性箱体的半尺寸（即半宽、半深和半高）为

\[\begin{align*} r_x = \sqrt{\frac{3}{2 \mathcal{M}} \left(\mathcal{I}_{yy} + \mathcal{I}_{zz} - \mathcal{I}_{xx} \right)} \\ r_y = \sqrt{\frac{3}{2 \mathcal{M}} \left(\mathcal{I}_{zz} + \mathcal{I}_{xx} - \mathcal{I}_{yy} \right)} \\ r_z = \sqrt{\frac{3}{2 \mathcal{M}} \left(\mathcal{I}_{xx} + \mathcal{I}_{yy} - \mathcal{I}_{zz} \right)} \end{align*} \]

设 \(\mathbf{v}\) 和 \(\boldsymbol{\omega}\) 分别表示物体在局部坐标系（与等效惯性箱体对齐）中的线速度和角速度。流体作用在固体上的力 \(\mathbf{f}_{\text{inertia}}\) 和扭矩 \(\mathbf{g}_{\text{inertia}}\) 是以下各项的总和

\[\begin{align*} \mathbf{f}_{\text{inertia}} &= \mathbf{f}_D + \mathbf{f}_V \\ \mathbf{g}_{\text{inertia}} &= \mathbf{g}_D + \mathbf{g}_V \end{align*} \]

此处的下标 \(D\) 和 \(V\) 分别表示二次阻力和粘性阻力。

二次阻力项取决于流体的密度 \(\rho\)，与物体速度的平方成比例，并且在高雷诺数下是液体力的有效近似。扭矩是通过对由旋转产生的力在表面积上积分而获得的。力和扭矩的第 \(i\) 个分量可以写成

\[\begin{aligned} f_{D, i} = \quad &- 2 \rho r_j r_k |v_i| v_i \\ g_{D, i} = \quad &- {1 \over 2} \rho r_i \left(r_j^4 + r_k^4 \right) |\omega_i| \omega_i \\ \end{aligned} \]

粘性阻力项取决于流体粘度 \(\beta\)，与物体速度呈线性关系，并近似低雷诺数下的液体力。注意，粘度可以独立于密度使用，以使仿真更阻尼。我们在低雷诺数下使用半径为 \(r_{eq} = (r_x + r_y + r_z) / 3\) 的等效球体的公式。局部物体坐标系中产生的 3D 力和扭矩为

\[\begin{aligned} f_{V, i} = \quad &- 6 \beta \pi r_{eq} v_i \\ g_{V, i} = \quad &- 8 \beta \pi r_{eq}^3 \omega_i \\ \end{aligned} \]

通过指定非零的风 (wind) 也可以影响这些力，风是在液体力学计算中从物体线速度中减去的 3D 向量。

椭球体模型#

../_images/fruitfly.png — 图中所示的可飞行果蝇 (Drosophila Melanogaster) 模型在Vaxenburg *et al.* [VSM+24]中有详细描述。#

在本节中，我们基于对几何形状的椭球体近似，描述并推导了一个无状态模型，用于模拟周围流体作用在运动刚体上的力。与上一节基于惯性的模型相比，该模型提供了对不同类型液体力的更精细控制。该模型的激励用例是昆虫飞行，参见右图。

概述#

通过将 fluidshape 属性设置为 ellipsoid 来为每个几何体激活该模型，这也会禁用父物体的基于惯性的模型。fluidcoef 属性中的 5 个数字对应于以下语义

索引	描述	符号	默认值
0	钝体阻力系数	\(C_{D, \text{blunt}}\)	0.5
1	细长体阻力系数	\(C_{D, \text{slender}}\)	0.25
2	角阻力系数	\(C_{D, \text{angular}}\)	1.5
3	库塔升力系数	\(C_K\)	1.0
4	马格努斯升力系数	\(C_M\)	1.0

该模型的元素是将 Andersen 等人 [APW05b] 推广到 3 维。流体作用在固体上的力 \(\mathbf{f}_{\text{ellipsoid}}\) 和扭矩 \(\mathbf{g}_{\text{ellipsoid}}\) 是以下各项的总和

\[\begin{align*} \mathbf{f}_{\text{ellipsoid}} &= \mathbf{f}_A + \mathbf{f}_D + \mathbf{f}_M + \mathbf{f}_K + \mathbf{f}_V \\ \mathbf{g}_{\text{ellipsoid}} &= \mathbf{g}_A + \mathbf{g}_D + \mathbf{g}_V \end{align*} \]

其中下标 \(A\)、\(D\)、\(M\)、\(K\) 和 \(V\) 分别表示附加质量、粘性阻力、马格努斯升力、库塔升力和粘性阻力。\(D\)、\(M\) 和 \(K\) 项分别按上述 \(C_D\)、\(C_M\) 和 \(C_K\) 系数缩放，粘性阻力随流体粘度 \(\beta\) 缩放，而附加质量项无法缩放。

符号#

我们描述物体在密度为 \(\rho\) 的无粘、不可压缩静止流体中的运动。模型中，任意形状的物体被描述为半轴为 \(\mathbf{r} = \{r_x, r_y, r_z\}\) 的等效椭球体。问题在一个与椭球体侧面平行并随之移动的参考系中描述。物体具有速度 \(\mathbf{v} = \{v_x, v_y, v_z\}\) 和角速度 \(\boldsymbol{\omega} = \{\omega_x, \omega_y, \omega_z\}\)。我们还将使用

\[\begin{align*} r_\text{max} &= \max(r_x, r_y, r_z) \\ r_\text{min} &= \min(r_x, r_y, r_z) \\ r_\text{mid} &= r_x + r_y + r_z - r_\text{max} - r_\text{min} \end{align*} \]

雷诺数是流体内惯性力与粘性力之比，定义为 \(Re=u~l/\beta\)，其中 \(\beta\) 是流体的运动粘度，\(u\) 是流动的特征速度（或者，通过坐标系变换，是物体的速度），\(l\) 是流动或物体的特征尺寸。

我们将使用 \(\Gamma\) 表示环量，它是速度场沿闭合曲线的线积分 \(\Gamma = \oint \mathbf{v} \cdot \textrm{d} \mathbf{l}\)，并且根据斯托克斯定理，\(\Gamma = \int_S \nabla \times \mathbf{v} \cdot \textrm{d}\mathbf{s}\)。在流体动力学符号中，符号 \(\boldsymbol{\omega}\) 通常用于表示涡量，定义为 \(\nabla \times \mathbf{v}\)，而不是角速度。对于刚体运动，涡量是角速度的两倍。

最后，我们使用下标 \(i, j, k\) 来表示对称应用于 \(x, y, z\) 的方程组。例如 \(a_i = b_j + b_k\) 是以下 3 个方程的简写

\[\begin{align*} a_x &= b_y + b_z \\ a_y &= b_x + b_z \\ a_z &= b_x + b_y \end{align*} \]

椭球体投影#

我们给出以下结果。

引理

给定一个半轴为 \((r_x, r_y, r_z)\) 且与坐标轴 \((x, y, z)\) 对齐的椭球体，以及一个单位向量 \(\mathbf{u} = (u_x, u_y, u_z)\)，该椭球体投影到与 \(\mathbf{u}\) 垂直的平面上的面积为

\[A^{\mathrm{proj}}_{\mathbf{u}} = \pi \sqrt{\frac{r_y^4 r_z^4 u_x^2 + r_z^4 r_x^4 u_y^2 + r_x^4 r_y^4 u_z^2}{r_y^2 r_z^2 u_x^2 + r_z^2 r_x^2 u_y^2 + r_x^2 r_y^2 u_z^2}} \]

展开查看推导

引理的推导

椭圆的面积

以原点为中心的任何椭圆都可以用二次形式 \(\mathbf{x}^T Q \mathbf{x} = 1\) 来描述，其中 \(Q\) 是一个实对称正定 2x2 矩阵，定义了椭圆的方向和半轴长度，\(\mathbf{x} = (x, y)\) 是椭圆上的点。椭圆的面积由下式给出

\[A = \frac{\pi}{\sqrt{\det Q}} . \]

椭球体截面

首先，我们计算以原点为中心的椭球体与通过原点且单位法向量为 \(\mathbf{n} = (n_x, n_y, n_z)\) 的平面 \(\Pi_{\mathbf{n}}\) 相交形成的椭圆面积。设 \((r_x, r_y, r_z)\) 是椭球体的半轴长度。不失一般性，假设椭球体的轴与坐标轴对齐就足够了。椭球体可以描述为 \(\mathbf{x}^T Q \mathbf{x} = 1\)，其中 \(Q = \textrm{diag}\mathopen{}\left( \left. 1 \middle/ r_x^2 \right., \left. 1 \middle/ r_y^2 \right., \left. 1 \middle/ r_z^2 \right. \right)\mathclose{}\)，\(\mathbf{x} = (x, y, z)\) 是椭球体上的点。

我们通过旋转平面 \(\Pi_{\mathbf{n}}\) 和椭球体，使得旋转后平面的法向量沿 \(z\) 轴方向。这样，通过将 \(z\) 坐标设置为零，我们就能获得所需的交集。将沿 \(z\) 轴方向的单位向量记为 \(\mathbf{\hat{z}}\)，我们有

\[\begin{align*} \mathbf{n} \times \mathbf{\hat{z}} &= \sin\theta \, \mathbf{m}, \\ \mathbf{n} \cdot \mathbf{\hat{z}} &= \cos\theta , \end{align*} \]

其中 \(\mathbf{m}\) 是定义旋转轴的单位向量，\(\theta\) 是旋转角。我们可以重新排列这些量，得到构成旋转四元数所需的量，即

\[\begin{align*} \cos\frac{\theta}{2} &= \sqrt{\frac{1+\cos\theta}{2}} &= \sqrt{\frac{1 + \mathbf{n} \cdot \mathbf{\hat{z}}}{2}}, \\ \sin\frac{\theta}{2}\,\mathbf{m} &= \frac{\mathbf{n} \times \mathbf{\hat{z}}}{2\cos\frac{\theta}{2}} &= \frac{\mathbf{n} \times \mathbf{\hat{z}}}{\sqrt{2 (1 + \mathbf{n} \cdot \mathbf{\hat{z}})}} . \end{align*} \]

因此，旋转四元数 \(q = q_r + q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k}\) 由下式给出

\[q_r = \sqrt{\frac{1 + n_z}{2}}, \quad q_x = \frac{n_y}{\sqrt{2 \left(1+n_z\right)}}, \quad q_y = \frac{-n_x}{\sqrt{2 \left(1+n_z\right)}}, \quad q_z = 0 . \]

由此，旋转矩阵由下式给出

\[\def\arraystretch{1.33} \begin{align*} R &= \begin{pmatrix} 1 - 2 q_y^2 - 2 q_z^2 & 2 \left(q_x q_y - q_r q_z\right) & 2 \left(q_x q_z + q_r q_y\right) \\ 2 \left(q_x q_y + q_r q_z\right) & 1 - 2 q_x^2 - 2 q_z^2 & 2 \left(q_y q_z - q_r q_x\right) \\ 2 \left(q_x q_z - q_r q_y\right) & 2 \left(q_y q_z + q_r q_x\right) & 1 - 2 q_x^2 - 2 q_y^2 \end{pmatrix} \\ &= \begin{pmatrix} 1 - \left. n_x^2 \middle/ \left( 1+n_z \right) \right. & \left. -n_x n_y \middle/ \left( 1+n_z \right) \right. & -n_x \\ \left. -n_x n_y \middle/ \left( 1+n_z \right) \right. & 1 - \left. n_y^2 \middle/ \left( 1+n_z \right) \right. & -n_y \\ n_x & n_y & 1 - \left. \left( n_x^2 + n_y^2 \right) \middle/ \left( \vphantom{n_x^2} 1+n_z \right) \right. \end{pmatrix}, \end{align*} \]

旋转后的椭球体通过变换后的二次形式描述

\[\mathbf{x}^T Q' \mathbf{x} = \mathbf{x}^T \left( R^T Q R \right) \mathbf{x} = 1 . \]

根据上面椭圆面积的公式，对于 \(z=0\) 处的椭圆面积，我们需要

\[\begin{align*} Q'_{xx} &= \frac{1}{r_x^2} R_{xx}^2 + \frac{1}{r_y^2} R_{yx}^2 + \frac{1}{r_z^2} R_{zx}^2 , \\ Q'_{yy} &= \frac{1}{r_x^2} R_{xy}^2 + \frac{1}{r_y^2} R_{yy}^2 + \frac{1}{r_z^2} R_{zy}^2 , \\ Q'_{xy} &= \frac{1}{r_x^2} R_{xx} R_{xy} + \frac{1}{r_y^2} R_{yx} R_{yy} + \frac{1}{r_z^2} R_{zx} R_{zy} , \end{align*} \]

所需面积由下式给出

\[A^{\cap}_{\mathbf{n}} = \frac{\pi}{\sqrt{\vphantom{Q'^2_{xy}} \det Q'}} = \frac{\pi}{\sqrt{Q'_{xx} Q'_{yy} - Q'^2_{xy}}} = \frac{\pi r_x r_y r_z}{\sqrt{r_x^2 n_x^2 + r_y^2 n_y^2 + r_z^2 n_z^2}}, \]

其中上标 \(\cap\) 表示面积对应于与 \(\Pi_{\mathbf{n}}\) 相交处的椭圆。

投影椭圆

设 \(\mathbf{u} = (u_x, u_y, u_z)\) 是某个单位向量（在我们的上下文中，它是流体冲击椭球体的速度方向），\(\Pi_{\mathbf{u}}\) 是与 \(\mathbf{u}\) 垂直的平面。一般来说，将椭球体 \(\mathcal{E}\) 投影到平面 \(\Pi_{\mathbf{u}}\) 上形成的椭圆（记作 \(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}}\)）与椭球体 \(\mathcal{E}\) 和平面 \(\Pi_{\mathbf{u}}\) 相交形成的椭圆（记作 \(\mathcal{E}^{\cap}_{\mathbf{u}}\)）是不同的。

\(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}}\) 的一个重要性质是，在 \(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}}\) 上的每一点， \(\mathbf{u}\) 都与椭球体 \(\mathcal{E}\) 相切。

我们可以将 \(\mathcal{E}\) 视为单位球 \(\mathcal{S}\) 在拉伸变换 \(T = \mathrm{diag}(r_x, r_y, r_z)\) 下的像。此外，如果 \(\mathbf{\tilde{u}}\) 是 \(\mathcal{S}\) 的切向量，则其像 \(\mathbf{u}=T\mathbf{\tilde{u}}=(r_x \tilde{u}_x, r_y \tilde{u}_y, r_z \tilde{u}_z)\) 是椭球体的切向量。因此，椭圆 \(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}}\) 是圆 \(\mathcal{C}^{\cap}_{\mathbf{\tilde{u}}}\) 在 \(T\) 下的像，该圆是 \(\mathcal{S}\) 与 \(\Pi_{\mathbf{\tilde{u}}}\) 相交处的圆（对于球体，\(\mathcal{C}^{\cap}\) 和 \(\mathcal{C}^{\mathrm{proj}}\) 是重合的）。

设 \(\mathbf{\tilde{v}}\) 和 \(\mathbf{\tilde{w}}\) 是平面 \(\Pi_{\mathbf{\tilde{u}}}\) 中的一对正交向量，则 \(\mathbf{\tilde{u}} = \mathbf{\tilde{v}} \times \mathbf{\tilde{w}}\)。它们在 \(T\) 下的像分别是 \(\mathbf{v} = (r_x \tilde{v}_x, r_y \tilde{v}_y, r_z \tilde {v}_z)\) 和 \(\mathbf{w} = (r_x \tilde{w}_x, r_y \tilde{w}_y, r_z \tilde {w}_z)\)，并且它们在 \(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}}\) 的平面中保持正交。因此，椭圆 \(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}}\) 的一个（非单位）法向量由下式给出

\[\mathbf{N} = \mathbf{v} \times \mathbf{w} = (r_y r_z \tilde{u}_x, r_z r_x \tilde{u}_y, r_x r_y \tilde{u}_z) = \left( \frac{r_y r_z}{r_x} u_x, \frac{r_z r_x}{r_y} u_y, \frac{r_x r_y}{r_z} u_z \right). \]

这表明 \(\mathcal{E}^{\mathrm{proj}}_{\mathbf{u}} = \mathcal{E}^{\cap}_{\mathbf{n}}\)，其中 \(\mathbf{n} = \mathbf{N} / \left\Vert\mathbf{N}\right\Vert\)。其面积由前一节导出的公式给出，得到了上述结果。

附加质量#

对于在流体中运动的物体，附加质量或虚拟质量衡量了由于物体的运动而移动的流体的惯性。它可以从势流理论中导出（也就是说，即使是无粘流体也存在附加质量）。

参照Lamb [Lam32] 的第五章，由于流体从静止开始运动而施加在运动物体上的力 \(\mathbf{f}_{V}\) 和扭矩 \(\mathbf{g}_{V}\) 可以写成

\[\begin{align*} \mathbf{f}_{A} &= - \frac{\textrm{d}}{\textrm{d} t} \nabla_{\mathbf{v}} \mathcal{T} + \nabla_{\mathbf{v}} \mathcal{T} \times \boldsymbol{\omega} \\ \mathbf{g}_{A} &= - \frac{\textrm{d}}{\textrm{d} t} \nabla_{\boldsymbol{\omega}} \mathcal{T} + \nabla_{\mathbf{v}} \mathcal{T} \times \mathbf{v} + \boldsymbol{\omega} \times \nabla_{\boldsymbol{\omega}} \mathcal{T} \end{align*} \]

其中 \(\mathcal{T}\) 是仅考虑流体的动能。这些力通常被称为附加质量力或虚拟质量力，因为它们是由于被加速物体移动或偏转流体的惯性所致。事实上，对于具有恒定线速度的物体，这些力会减小到零。我们将物体视为具有三个对称平面，因为在此假设下，动能大大简化，可以写成

\[2 \mathcal{T} = m_{A, x} v_x^2 + m_{A, y} v_y^2 + m_{A, z} v_z^2 + I_{A, x} \omega_x^2 + I_ {A, y} \omega_y^2 + I_{A, y} \omega_z^2 \]

为了方便，我们引入附加质量向量 \(\mathbf{m}_A = \{m_{A, x}, m_{A, y}, m_{A, z}\}\) 和附加转动惯量向量 \(\mathbf{I}_A = \{I_{A, x}, I_{A, y}, I_{A, z}\}\)。这些量中的每一个都应该估计由于物体沿相应方向的运动而移动的流体的惯性，并且可以从势流理论中为一些简单的几何形状推导出来。

对于具有三个对称平面的物体，我们可以用紧凑的形式写出由附加惯性产生的力和扭矩

\[\begin{align*} \mathbf{f}_{A} &= - \mathbf{m}_A \circ \dot{\mathbf{v}} + \left(\mathbf{m}_A \circ \mathbf{v} \right) \times \boldsymbol{\omega} \\ \mathbf{g}_{A} &= - \mathbf{I}_A \circ \dot{\boldsymbol{\omega}} + \left(\mathbf{m}_A \circ \mathbf{v} \right) \times \mathbf{v} + \left(\mathbf{I}_A \circ \boldsymbol{\omega} \right) \times \boldsymbol{\omega} \end{align*} \]

这里 \(\circ\) 表示按元素乘积， \(\dot{\mathbf{v}}\) 是线加速度， \(\dot{\boldsymbol{\omega}}\) 是角加速度。 \(\mathbf{m}_A \circ \mathbf{v}\) 和 \(\mathbf{I}_A \circ \boldsymbol{\omega}\) 分别是虚拟线动量和虚拟角动量。

对于半轴为 \(\mathbf{r} = \{r_x, r_y, r_z\}\) 体积为 \(V = 4 \pi r_x r_y r_z / 3\) 的椭球体，虚拟惯性系数由Tuckerman [Tuc25] 推导。令

\[\kappa_i = \int_0^\infty \frac{r_i r_j r_k}{\sqrt{(r_i^2 + \lambda)^3 (r_j^2 + \lambda) (r_k^2 + \lambda)}} \textrm{d} \lambda \]

应该注意，这些系数是无量纲的（即如果所有半轴都乘以同一个标量，系数保持不变）。椭球体的虚拟质量是

\[m_{A, i} = \rho V \frac{\kappa_i}{2 - \kappa_i} \]

虚拟转动惯量是

\[I_{A, i} = \frac{\rho V}{5} \frac{(r_j^2 - r_k^2)^2 (\kappa_k-\kappa_j)}{2(r_j^2 - r_k^2) + (r_j^2 + r_k^2) (\kappa_j-\kappa_k)} \]

粘性阻力#

阻力作用于抵抗物体相对于周围流体的运动。我们发现粘性力也有助于降低包含流体动力学项的运动方程的刚度。因此，我们选择偏向保守的一面，并选择了可能高估耗散的粘性项近似值。

尽管最终是由粘性耗散引起，但在高雷诺数下，阻力与粘度无关，并与速度的平方成比例。它可以写成

\[\begin{align*} \mathbf{f}_\text{D} = - C_D~\rho~ A_D ~ \|\mathbf{v}\|~ \mathbf{v}\\ \mathbf{g}_\text{D} = - C_D \rho~ I_D ~ \|\boldsymbol{\omega}\| ~ \boldsymbol{\omega} \end{align*} \]

其中 \(C_D\) 是阻力系数， \(A_D\) 是参考表面积（例如，衡量垂直于流动的平面上的投影面积）， \(I_D\) 是参考转动惯量。

即使对于简单的形状，术语 \(C_D\)、\(A_D\) 和 \(I_D\) 也需要根据具体问题的物理和动力学尺度进行调整 [DHD15]。例如，阻力系数 \(C_D\) 通常随雷诺数的增加而减小，并且单个参考面积 \(A_D\) 可能不足以解释高度不规则或细长体的表面阻力。例如，实验拟合来自从下落的扑克牌 [APW05a, APW05b, WBD04] 到粒子输运 [BB16, Lot08] 的问题。参见右侧 cards.xml 模型的屏幕截图。

我们基于两个表面 \(A^\text{proj}_\mathbf{v}\) 和 \(A_\text{max}\) 推导了 \(\mathbf{f}_\text{D}\) 的公式。第一个表面 \(A^\text{proj}_\mathbf{v}\) 是物体在垂直于速度 \(\mathbf{v}\) 的平面上的圆柱投影。第二个表面是最大投影面积 \(A_\text{max} = 4 \pi r_{max} r_{min}\)。

\[\mathbf{f}_\text{D} = - \rho~ \big[ C_{D, \text{blunt}} ~ A^\text{proj}_\mathbf{v} ~ + C_{D, \text{slender}}\left(A_\text{max} - A^\text{proj}_\mathbf{v} \right) \big] ~ \|\mathbf{v}\|~ \mathbf{v} \]

\(A^\text{proj}_\mathbf{v}\) 的公式和推导在上面的引理中给出。

我们为角阻力提出了一个类似的模型。对于每个笛卡尔轴，我们考虑物体围绕该轴旋转时获得的最大扫掠椭球体的转动惯量。由此产生的转动惯量对角项是

\[\mathbf{I}_{D,ii} = \frac{8\pi}{15} ~r_i ~\max(r_j, ~r_k)^4 . \]

给定这个参考转动惯量，角阻力扭矩计算为

\[\mathbf{g}_\text{D} = - \rho ~ \boldsymbol{\omega} ~ \Big( \big[ C_{D, \text{angular}} ~ \mathbf{I}_D ~ + C_{D, \text{slender}} \left(\mathbf{I}_\text{max} - \mathbf{I}_D \right) \big] \cdot \boldsymbol{\omega} \Big) \]

这里 \(\mathbf{I}_\text{max}\) 是一个向量，其每个分量等于 \(\mathbf{I}_D\) 的最大分量。

最后，粘性阻力项，也称为线性阻力，在雷诺数约为或低于 \(O(10)\) 时能很好地近似流体作用力。这些力是根据等效球体的斯托克斯定律计算的 [Lam32, Sto50]

\[\begin{align*} \mathbf{f}_\text{V} &= - 6 \pi r_D \beta \mathbf{v}\\ \mathbf{g}_\text{V} &= - 8 \pi r_D^3 \beta \boldsymbol{\omega} \end{align*} \]

这里， \(r_D = (r_x + r_y + r_z)/3\) 是等效球体的半径， \(\beta\) 是介质的运动粘度（例如，常温空气约为 \(1.48~\times 10^{-5}~m^2/s\)，水约为 \(0.89 \times 10^{-4}~m^2/s\)）。举一个定量的例子，如果 \(u\cdot l \lesssim 2 \times 10^{-4}~m^2/s\)，斯托克斯定律对于室温空气是精确的，其中 \(u\) 是速度， \(l\) 是物体的特征长度。

粘性升力#

库塔-儒可夫斯基定理计算了一个在均匀流中以速度 u 平移的二维物体的升力 \(L\) 为 \(L = \rho u \Gamma\)。这里， \(\Gamma\) 是物体周围的环量。在接下来的小节中，我们定义了两种环量的来源以及由此产生的升力。

马格努斯力#

../_images/magnus.png — 绕旋转圆柱体的流体烟流可视化 (WikiMedia Commons, CC BY-SA 4.0)。由于粘度，旋转的圆柱体使来流向上偏转，并受到向下的力（红色箭头）。#

马格努斯效应描述了旋转物体在流体中运动时的现象。通过粘性效应，旋转物体会在周围流体中诱导旋转。这种旋转会使流体绕过物体的轨迹发生偏转（即导致线加速度），而物体会受到一个相等且相反的反作用力。对于圆柱体，单位长度圆柱体的马格努斯力可以计算为 \(F_\text{M} / L = \rho v \Gamma\)，其中 \(\Gamma\) 是由旋转引起的流体环量， \(v\) 是物体的速度。我们对于任意物体估计这个力为

\[\mathbf{f}_{\text{M}} = C_M ~\rho~ V~ \boldsymbol{\omega}\times\mathbf{v} , \]

其中 \(V\) 是物体的体积， \(C_M\) 是力的系数，通常设为 1。

举个例子很有意义。为了减少变量，假设物体只沿一个方向旋转，例如 \(\boldsymbol{\omega} = \{0, 0, \omega_z\}\)，并沿另外两个方向平移，例如 \(\mathbf{v} = \{v_x, v_y, 0\}\)。附加质量力与马格努斯效应力沿，例如， \(x\) 方向的合力是

\[\frac{f}{\pi \rho r_z} = v_y \omega_z \left(2 r_x \min\{r_x, r_z\} - (r_x + r_z)^2\right) \]

注意这两个项的符号相反。

库塔条件#

驻点是流场中速度为零的位置。对于在流体中运动的物体（在二维情况下，在随物体运动的参考系中），存在两个驻点：一个在前端，流线在此分离并流向物体的两侧；另一个在后端，流线在此重新汇合。一个带有尖锐尾缘（后端边缘）的运动物体会在周围流体中产生足够强的环量，以将后驻点保持在尾缘。这就是库塔条件，一个流体动力学现象，可以在带有尖角的固体物体上观察到，例如细长体或翼型的后缘。

../_images/kutta_cond_plate_dark.svg — 库塔条件示意图。蓝色线条是流线，两个品红色点是驻点。连接两个驻点的分离流线用绿色标记。分离流线和物体划定了一个区域，据说该区域流体“分离”并在内部再循环。这种环量在平板上产生了向上的力。#

对于上面图中所示的二维流场，由于库塔条件产生的环量可以估计为： \(\Gamma_\text{K} = C_K ~ r_x ~ \| \mathbf{v}\| ~ \sin(2\alpha)\)，其中 \(C_K\) 是升力系数， \(\alpha\) 是速度向量与其在表面上投影之间的夹角。单位长度的升力可以根据库塔-儒可夫斯基定理计算为 \(\mathbf{f}_K / L = \rho \Gamma_\text{K} \times \mathbf{v}\)。

为了将升力方程扩展到三维运动，我们考虑物体横截面的法向量 \(\mathbf{n}_{s, \mathbf{v}} = \{\frac{r_y r_z}{r_x}v_x, \frac{r_z r_x}{r_y}v_y, \frac{r_x r_x}{r_z}v_z\}\)，该横截面产生了物体在垂直于速度的平面上的投影 \(A^\text{proj}_\mathbf{v}\)（如上面的引理中所示），以及相应的单位向量 \(\hat{\mathbf{n}}_{s, \mathbf{v}}\)。我们使用这个方向将 \(\mathbf{v}\) 分解为 \(\mathbf{v}_\parallel ~+~ \mathbf{v}_\perp\)，其中 \(\mathbf{v}_\perp = \left(\mathbf{v} \cdot \hat{\mathbf{n}}_{s, \mathbf{v}}\right) \hat{\mathbf{n}}_{s, \mathbf{v}}\)。我们写出升力为

\[\begin{align*} \mathbf{f}_\text{K} &= \frac{C_K~\rho~ A^\text{proj}_\mathbf{v}}{\|\mathbf{v}\|} \left( \mathbf{v} \times \mathbf{v}_\parallel\right)\times \mathbf{v} \\ &= C_K~\rho~ A^\text{proj}_\mathbf{v} \left(\hat{\mathbf{v}} \cdot \hat{\mathbf{n}}_{s, \mathbf{v}}\right) \left( \hat{\mathbf{n}}_{s, \mathbf{v}} \times \mathbf{v} \right)\times \mathbf{v} \end{align*} \]

这里， \(\hat{\mathbf{v}}\) 是沿 \(\mathbf{v}\) 方向的单位法向量。注意， \(\hat{\mathbf{n}}_{s, \mathbf{v}}\) 的方向仅在物体半轴不相等的平面上与 \(\hat{\mathbf{v}}\) 不同。因此，例如，对于球形物体， \(\hat{\mathbf{n}}_{s, \mathbf{v}} \equiv \hat{\mathbf{v}}\)，并且根据定义 \(\mathbf{f}_\text{K} = 0\)。

让我们用一个例子来阐述这个关系。假设一个物体具有 \(r_x = r_y\) 和 \(r_z \ll r_x\)。注意向量 \(\hat{\mathbf{n}}_{s, \mathbf{v}} \times \hat{\mathbf{v}}\) 给出了由固体物体偏转流体引起的环量方向。沿 \(z\) 方向，环量将正比于 \(\frac{r_y r_z}{r_x}v_x v_y - \frac{r_z r_x}{r_y}v_x v_y = 0\)（由于 \(r_x = r_y\)）。因此，在固体形状钝圆的平面上，运动不会产生环量。

现在，为了简化，令 \(v_x = 0\)。在这种情况下，沿 \(y\) 方向的环量也为零，它正比于 \(\frac{r_y r_z}{r_x}v_x v_z - \frac{r_y r_x}{r_y}v_x v_z\)。唯一的非零环量分量将沿 \(x\) 方向，并正比于 \(\left(\frac{r_x r_z}{r_y} - \frac{r_x r_y}{r_z}\right) v_y v_z \approx \frac{r_x^2}{r_z} v_y v_z\)。

我们将有 \(\mathbf{v}_\parallel = \{v_x, 0, v_z\}\)，环量 \(\Gamma \propto \{r_z v_y v_z, ~ 0,~ - r_x v_x v_y \} / \|\mathbf{v}\|\)。运动在固体形状钝圆的平面上不会产生环量，而在另外两个平面上，环量为 \(\Gamma \propto r_\Gamma ~ \|\mathbf{v}\|~ \sin(2 \alpha) ~ = ~2 r_\Gamma ~\|\mathbf{v}\| ~\sin(\alpha)~\cos(\alpha)\)，其中 \(\alpha\) 是速度与其在平面上（例如，垂直于 \(x\) 的平面上）物体投影之间的夹角（例如，在垂直于 \(x\) 的平面上，我们有 \(\sin(\alpha) = v_y/\|\mathbf{v}\|\) 和 \(\cos(\alpha) = v_z/\|\mathbf{v}\|\)），而 \(r_\Gamma\) 是平面上的升力表面（例如，对于垂直于 \(x\) 的平面，是 \(r_z\)）。此外，环量的方向由叉积给出（因为固体边界将输入流速“旋转”到其在物体上的投影方向）。

致谢#

本节中模型的设计和实现是 Guido Novati 的工作。

参考文献#

APW05a: A Andersen, U Pesavento, and Z Jane Wang. Unsteady aerodynamics of fluttering and tumbling plates. Journal of Fluid Mechanics, 541:65–90, 2005.
APW05b(1,2): Anders Andersen, Umberto Pesavento, and Z Jane Wang. Analysis of transitions between fluttering, tumbling and steady descent of falling cards. Journal of Fluid Mechanics, 541:91–104, 2005.
BB16: Gholamhossein Bagheri and Costanza Bonadonna. On the drag of freely falling non-spherical particles. Powder Technology, 301:526–544, 2016.
DHD15: Zhipeng Duan, Boshu He, and Yuanyuan Duan. Sphere drag and heat transfer. Scientific reports, 5(1):1–7, 2015.
Lam32(1,2): Horace Lamb. Hydrodynamics. Sixth edition. Cambridge University Press, 1932.
Lot08: E Loth. Drag of non-spherical solid particles of regular and irregular shape. Powder Technology, 182(3):342–353, 2008.
Sto50: GG Stokes. On the effect of internal friction of fluids on the motion of pendulums. Trans. Camb. phi1. S0c, 9(8):106, 1850.
Tuc25: LB Tuckerman. Inertia factors of ellipsoids for use in airship design. US Government Printing Office, 1925.
VSM+24: Roman Vaxenburg, Igor Siwanowicz, Josh Merel, Alice A Robie, Carmen Morrow, Guido Novati, Zinovia Stefanidi, Gwyneth M Card, Michael B Reiser, Matthew M Botvinick, Kristin M Branson, Yuval Tassa, and Srinivas C Turaga. Whole-body simulation of realistic fruit fly locomotion with deep reinforcement learning. bioRxiv, 2024. doi:10.1101/2024.03.11.584515.
WBD04: Z Jane Wang, James M Birch, and Michael H Dickinson. Unsteady forces and flows in low reynolds number hovering flight: two-dimensional computations vs robotic wing experiments. Journal of Experimental Biology, 207(3):449–460, 2004.