Cosserat Rod 理论学习

前言

最近做的本科毕设中有一部分的工作是要对缝合线进行仿真，老师推荐用 Cosserat Rod 模型对缝合线建模，于是我就查阅了很多关于 Cosserat Rod 的资料，并以本文作为学习笔记。

本文首先介绍了描述空间曲线的 TNB Frame 和 Darboux Vector，然后以此为基础介绍了 Cosserat Rod 模型的相关理论，最后根据 Position and Orientation Based Cosserat Rods 这篇论文介绍了 Cosserat Rod 的离散形式以及如何用 PBD 算法进行仿真。

TNB Frame

大多数的对缝合线、毛发、绳子等物体仿真使用的模型都是线型模型，也就是要对空间曲线的运动进行模拟。因此首先介绍一下 TNB Frame 的相关理论。

TNB Frame（中文翻译应该是 TNB 框架或者 TNB 参考系？）也叫作 Frenet–Serret frame，是由这两个人分别独立研究提出来的理论，所以以这两个人的名字命名。TNB Frame 理论描述了一个粒子沿着三维空间中连续可微曲线的运动。所谓的 “TNB”，指的是三维空间中的一组正交基，即三个相互垂直的单位向量：

• T：Tangent，曲线的单位切向量；
• N：Normal，曲线的单位法向量；
• B：Binormal，曲线的单位副法向量；

假设空间曲线的方程定义为 $\mathbf{r}(s)$，其中 $s$ 表示弧长且 $s\in [0,L]$，那么 TNB 三个向量可以通过 $\mathbf{r}(s)$ 计算得到：

首先单位切向量 $\mathbf{T}=\frac{d\mathbf{r}}{ds}/||\frac{d\mathbf{r}}{ds}||$ ，它指向了粒子在曲线当前位置运动的前进方向，而 $\frac{d\mathbf{T}}{ds}$ 则可以表示这个粒子的加速度的方向，注意我们这里的速度和加速度的定义不是严格的物理意义的定义，而是为了方便描述曲线才假设有个粒子在曲线上面运动，并给粒子赋予了速度和加速度，用同样为标量的弧长 $s$ 的变化替代时间变化。曲线的曲率 $\kappa (s)=||\frac{d\mathbf{T}}{ds}||$，则单位法向量 $\mathbf{N}=\frac{1}{\kappa }\frac{d\mathbf{T}}{ds}$ 。TNB 三个向量组成了正交基，因此单位副法向量 $\mathbf{B}=\mathbf{T}\times \mathbf{N}$ 。总体来说，假设有一个粒子沿着曲线匀速运动，那么 $\mathbf{T}$ 就是这个粒子的速度方向，$\mathbf{N}$ 是这个粒子的加速度方向，$\mathbf{B}$ 则垂直于上面两个方向。总结一下：

• $\mathbf{T}=\frac{d\mathbf{r}}{ds}/||\frac{d\mathbf{r}}{ds}||$
• $\mathbf{N}=\frac{1}{\kappa }\frac{d\mathbf{T}}{ds}$，其中曲率 $\kappa (s)=||\frac{d\mathbf{T}}{ds}||$
• $\mathbf{B}=\mathbf{T}\times \mathbf{N}$

TNB Frame 的相关理论中有一套公式 Frenet–Serret formulas：
$$\begin{array}{rl}\frac{d\mathbf{T}}{ds}& =\kappa \mathbf{N}\\ \frac{d\mathbf{N}}{ds}& =-\kappa \mathbf{T}+\tau \mathbf{B}\\ \frac{d\mathbf{B}}{ds}& =-\tau \mathbf{N}\end{array}$$
其中 $\kappa$ 在前面说过了表示的是曲率（curvature），$\tau$ 表示的是曲线的扭率（torsion），两者都可以反推回来计算出来：
$$\begin{array}{rl}\kappa (s)& =||\frac{d\mathbf{T}}{ds}||\\ \tau (s)& =-\mathbf{N}\cdot \frac{d\mathbf{B}}{ds}\end{array}$$
其中 $\tau$ 的计算公式成立是因为 $\mathbf{N}$ 为单位向量。从直观上来说，曲率 $\kappa$ 刻画了曲线的弯曲程度（刻画了 $\mathbf{B}$ 的转动），扭率 $\tau$ 刻画了曲线的扭曲程度（刻画了 $\mathbf{T}$ 的转动）。

上面的叙述可能有点乱，但是从实际运用的角度来说，一般情况下空间曲线 $\mathbf{r}(s)$ 是已知的，其他的相关量 { T , N , B , κ , τ } \{\boldsymbol{T},\boldsymbol{N},\boldsymbol{B},\kappa,\tau\} {T,N,B,κ,τ} 都可以计算出来。

Darboux Vector

达布向量 $\mathbf{\omega }$ 描述的是 TNB Frame 中三个单位向量基的转动，定义为 $\mathbf{\omega }={\mathbf{\omega }}_{\mathbf{T}}+{\mathbf{\omega }}_{\mathbf{N}}+{\mathbf{\omega }}_{\mathbf{B}}$ ，其中后面三个 ${\mathbf{\omega }}_{\mathbf{x}}$ 分别描述 TNB 三个单位向量的略面速度，例如
$${\mathbf{\omega }}_{\mathbf{T}}=\underset{\Delta t\to 0}{\lim }\frac{\mathbf{T}(t)\times \mathbf{T}(t+\Delta t)}{2\Delta t}=\frac{\mathbf{T}\times {\mathbf{T}}^{\prime }}{2}$$
根据定义可以将 $\mathbf{\omega }$ 的计算公式简化得到
$$\mathbf{\omega }=\tau \mathbf{T}+\kappa \mathbf{B}$$
达布向量具有如下性质：
$$\begin{array}{rl}{\mathbf{T}}^{\prime }& =\mathbf{\omega }\times \mathbf{T}\\ {\mathbf{N}}^{\prime }& =\mathbf{\omega }\times \mathbf{N}\\ {\mathbf{B}}^{\prime }& =\mathbf{\omega }\times \mathbf{B}\end{array}$$
结合我们熟知的公式 $\mathbf{v}=\mathbf{\omega }\times \mathbf{r}$ 来看达布向量可以发现，达布向量类似于 TNB Frame 里面的角速度。

Cosserat Rod

Cosserat Rod（图片来自文献Position and Orientation Based Cosserat Rods）

在 Cosserat Rod 理论中认为线型模型的长度远远大于横截面半径，即 $L>>r$，在使用时直接认为是一段空间曲线。有两个坐标系，一个是全局的坐标系 { e 1 , e 2 , e 3 } \{\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3\} {e1​,e2​,e3​}，另一个是与弧长 $s$ 相关的局部坐标系 { d 1 ( s ) , d 2 ( s ) , d 3 ( s ) } \{\boldsymbol{d}_1(s),\boldsymbol{d}_2(s),\boldsymbol{d}_3(s)\} {d1​(s),d2​(s),d3​(s)} ，在局部坐标系中， { d 1 ( s ) , d 2 ( s ) } \{\boldsymbol{d}_1(s),\boldsymbol{d}_2(s)\} {d1​(s),d2​(s)} 构成了曲线上某一个点的横截面，而 ${\mathbf{d}}_3(s)$ 是这个横截面的法向量，注意这个 ${\mathbf{d}}_3(s)$ 不一定要和曲线切线 $\frac{d\mathbf{r}(s)}{ds}$ 同方向，当它们方向不同时说明在 $s$ 处有剪切。

可以看出 Cosserat Rod 局部坐标系的定义和 TNB Frame 很相似，唯一不同的就是它不严格要求空间曲线一定连续可微。

关于全局坐标系和局部坐标系的变换，我们可以看成是一种旋转，定义表示旋转的四元数 $q(s)$，则有
$$d_k=q{e}_k{q}^{\ast }$$
Cosserat Rod 模型的好处在于，它可以度量线型模型的延展（stretch）、剪切（shear）、弯曲（bend）和扭转（twist）四种应变，因此可以更加真实地模拟线型物体的物理运动。

定义
$$\mathbf{\Gamma }(s)=\frac{d\mathbf{r}(s)}{ds}-{\mathbf{d}}_3(s)$$
用于度量曲线的延展和剪切应变，${\mathbf{d}}_3(s)$ 和 $\frac{d\mathbf{r}(s)}{ds}$ 的方向关系表示了其剪切应变，${\mathbf{d}}_3(s)$ 和 $\frac{d\mathbf{r}(s)}{ds}$ 的长度关系表示了其延展应变。

定义
$$\Delta \omega =\omega -{\omega }^0$$
用于度量曲线的弯曲和扭转应变，其中 $\omega$ 对应达布向量 $\mathbf{\omega }$ 的四元数，上标 $0$ 表示初始值。在上一个部分中我们知道了达布向量的计算方法
$$\mathbf{\omega }=\frac{1}{2}\sum _{k=1}^3{\mathbf{d}}_k\times {\mathbf{d}}_k^{\prime }$$
也知道了达布向量 $\mathbf{\omega }$ 描述了 { d 1 ( s ) , d 2 ( s ) , d 3 ( s ) } \{\boldsymbol{d}_1(s),\boldsymbol{d}_2(s),\boldsymbol{d}_3(s)\} {d1​(s),d2​(s),d3​(s)} 随着 $s$ 的变化的角速度，因此可以发现达布向量的前两个元素 { ω 1 , ω 2 } \{\omega_1,\omega_2\} {ω1​,ω2​} 描述了曲线的弯曲应变（影响的是曲率），而最后一个元素 ${\omega }_3$ 描述了曲线的扭转应变（影响的是扭率）。

结合局部坐标系和全局坐标系的变换四元数 $q$ ，我们可以将两个应变度量函数改写为
$$\begin{array}{rl}\mathbf{\Gamma }(s)& =\mathbf{R}(q{)}^T\frac{d\mathbf{r}(s)}{ds}-{\mathbf{e}}_3\\ \Delta \omega & =2(q^{\ast }{q}^{\prime }-q^{\ast 0}{q}^{\prime 0})\end{array}\text{\hspace{1em}(1)}$$
其中 $\mathbf{R}(q)$ 为四元数 $q$ 对应的旋转矩阵（$qv{q}^{\ast }$ 与 $\mathbf{R}(q)\mathbf{v}$ 表示同一种旋转变换，相互可以推导，过程比较复杂所以这里省略了），$q^{\prime }$ 为四元数的导数，有 $q^{\prime }=\frac{1}{2}q\omega$ 。

离散形式的 Cosserat Rod

物理仿真算法中都是对离散节点的运动进行仿真，因此 Position and Orientation Based Cosserat Rods 这篇论文构建了离散形式的 Cosserat Rod，并且成功套用了 PBD 的约束作用于节点位置和节点间的四元数。

离散形式的 Cosserat Rod（图片来自文献Position and Orientation Based Cosserat Rods）

假设曲线由 n 个节点 ${\mathbf{x}}_k$ 组成，那么相邻节点之间还添加了一个表示一段杆（rod）的节点 $q_{k-\frac{1}{2}}$ 。

在杆节点处构建约束最小化延展和剪切应变，在位置节点处构建约束最小化弯曲和扭转，于是我们需要把原本公式 (1) 中的一些量进行离散化处理。对于 $\frac{d\mathbf{r}(s)}{ds}$ 离散化为 $\frac{d\mathbf{r}(s)}{ds}\approx \frac{1}{l}({\mathbf{x}}_{i+1}-{\mathbf{x}}_i)$ ，其中 $l$ 表示杆节点的长度；对于 $\mathbf{\omega }$ 离散化为 $\mathbf{\omega }\approx \frac{2}{l}\varsigma (q_{i-\frac{1}{2}}^{\ast }{q}_{i+\frac{1}{2}})$ ，其中 $\varsigma$ 表示取四元数的向量部分，这个离散化的推导我看了相关论文还不明白，之后有时间再慢慢去学习吧。

离散化之后的方程变为
$$\begin{array}{rl}{\mathbf{\Gamma }}_{i+\frac{1}{2}}& \approx \frac{1}{l}({\mathbf{x}}_{i+1}-{\mathbf{x}}_i)-\varsigma (q_{i+\frac{1}{2}}{e}_3{q}_{i+\frac{1}{2}}^{\ast })\\ \Delta \mathbf{\omega }& \approx \frac{2}{l}\varsigma (q_{i-\frac{1}{2}}^{\ast }{q}_{i+\frac{1}{2}}-q_{i-\frac{1}{2}}^{\ast 0}{q}_{i+\frac{1}{2}}^0)\end{array}\text{\hspace{1em}(2)}$$
那么对于 PBD 中如何构建约束，只要把位置 ${\mathbf{x}}_i$ 和表示转动的四元数 $q_{i-\frac{1}{2}}$ 当成被约束量，根据公式 (2) 令它们等于 $0$ 即可：
$$\begin{array}{rl}{\mathbf{C}}_{\mathbf{s}}({\mathbf{p}}_1,{\mathbf{p}}_2,q)& =\frac{1}{l}({\mathbf{p}}_2-{\mathbf{p}}_1)-\mathbf{R}(q){\mathbf{e}}_3=0\\ {\mathbf{C}}_{\mathbf{b}}(q,u)& =\varsigma (q^{\ast }u-q^{\ast 0}{u}^0)=\mathbf{\omega }-s{\mathbf{\omega }}^0\end{array}\text{\hspace{1em}(3)}$$
其中
$$s=\{\begin{array}{lr}+1& \text{if }|\mathbf{\omega }-{\mathbf{\omega }}^0{|}^2<|\mathbf{\omega }+{\mathbf{\omega }}^0{|}^2\\ -1& \text{if }|\mathbf{\omega }-{\mathbf{\omega }}^0{|}^2>|\mathbf{\omega }+{\mathbf{\omega }}^0{|}^2\end{array}$$
因为四元数 $+q$ 和 $-q$ 描述的是同一种转动，所以要加一个 $s$ 变量来确定初始值的唯一性。

根据 PBD 算法的公式可以得到每一帧的更新方程，对于约束 ${\mathbf{C}}_{\mathbf{s}}$，更新方程为
$$\begin{array}{rl}\Delta {\mathbf{p}}_1& =+\frac{w_{1}l}{w_1+w_2+4w_q{l}^2}(\frac{1}{l}({\mathbf{p}}_2-{\mathbf{p}}_1)-{\mathbf{d}}_3)\\ \Delta {\mathbf{p}}_2& =-\frac{w_{2}l}{w_1+w_2+4w_q{l}^2}(\frac{1}{l}({\mathbf{p}}_2-{\mathbf{p}}_1)-{\mathbf{d}}_3)\\ \Delta q& =+\frac{w_{q}l}{w_1+w_2+4w_q{l}^2}(\frac{1}{l}({\mathbf{p}}_2-{\mathbf{p}}_1)-{\mathbf{d}}_3)q{e}_3^{\ast }\end{array}$$
对于约束 ${\mathbf{C}}_{\mathbf{b}}$ ，更新方程为
$$\begin{array}{rl}\Delta q& =+\frac{w_q}{w_q+w_u}u(\omega -s{\omega }^0)\\ \Delta u& =-\frac{w_u}{w_q+w_u}q(\omega -s{\omega }^0)\end{array}$$