横向各向同性软组织的本构方程及有限元实现

前言

软组织的力学研究在近几十年内已经成熟，各种理论层出不穷，但总体来讲经历了三个阶段：

第一阶段，1940年Mooney提出了一个具有很强变形能力的超弹性材料的力学模型，1948年Rivlin对它做了改良，也就是现在流行的Mooney-Rivlin超弹性模型，目前对于超弹性材料的仿真都使用了Mooney-Rivlin。

第二阶段，1994年，Jeffrey Weiss 定义了纤维增强型Mooney-Rivlin，通过向量场控制当前构型的变形，这个理论为肌肉的力学仿真奠定了基础，并发表了论文（他博士读了九年，都在搞这个）。

第三阶段，2003年-2018年，在斯坦福大学跟工业光魔的共同努力下，肌肉仿真成功在数字人上应用，大幅提升了数字形象的真实感，目前 Ziva Dynamics，FACEGOOD 等公司也是类似技术路线。

我在2016年看到了Weta Digital一篇专访，放了这段肌肉仿真的视频（ Weta Digital Horse FEM Simulation）它们从内到外做了一匹“真马”：骨骼层->肌肉层->筋膜层->脂肪层->皮层，深深被其震撼，了解到他们用到了有限元法（FEM），随后就沉迷这个技术了。我个人来讲喜欢写一些有难度的玩意，编译器、虚拟机、3D引擎、渲染器、物理引擎之类的都写过，这么牛逼的物理引擎不搞一下自然心痒得很，很幸运，是兴趣又跟自己的创业方向吻合，肌肉仿真的研发就这样开始了。刚开始对计算力学一窍不通，只好从头学起，一开始就像个无头苍蝇学习路线比较乱，有限元的书买了一堆，发现它只是个数值工具，想搞明白原理需要先看力学，力学里面又有经典力学、理论力学、固体力学、结构力学、材料力学、弹性力学，不知道从哪里出发，于是干脆全都看了一遍，最后范围逐渐缩小，聚焦之后学习起来自然就很高效了，现在已经在手撸有限元的C++实现了。

图（2）计算力学知名学者吴建营老师（左）

2017年FACEGOOD上线了表情捕捉产品，但当时对计算力学方向还是很迷惘，有幸遇见了吴建营老师，我把这个课题给他看了之后他给的反馈是“都是技术问题，没有理论问题，但是工程化方面人才比较缺”。听了这一番话我反倒是觉着靠谱了，工程化是我们的专长，这些年一直在做工程化的急先锋。FACEGOOD 搞仿真主要是提升数字人的真实感，例如只有从力学上模拟表情产生的全过程才能做到Physically Correct，这里面最核心的是解决人体软组织的力学建模问题，也就是需要得到韧带、肌腱、脂肪、皮层等软组织的精确力学模型。本文主要介绍软组织的数学建模以及其有限元的实现，读者需要具备变分法（Caculus of Variations）、计算力学（Computional Mechanics）、有限元法（Finite Element Method）等知识，熟悉横向各向同性（Transverse Isotropy）、不可压缩超弹性材料的应力、弹性张量等理论特性。按照本栏目一向硬核风格，直接给出完整思路，最后附上一段小DEMO。

认识韧带与肌腱

图（3）

人体有206块骨头，大部分骨头是不能直接向接触的，这样会有很大的摩擦，也没法运动。骨与骨之间并没有榫卯结构，也没有机械齿轮，把骨与骨连结起来的就是韧带和关节。从解剖学角度讲，韧带是骨与骨之间的直接连结，指的就是连接骨与骨之间的致密结缔组织。

每一块骨骼肌都分成肌腹和肌腱两部分，肌腹由肌纤维构成，色红质软，有收缩能力，肌腱由致密结缔组织构成，色白较硬，没有收缩能力。肌腱为肌肉末端的结缔组织纤维索，肌肉藉此附着于骨骼或其它结构。

材料初步（超弹性，横向各向异性的Strain Energy方程推导）

大部分生物材料的力学特性都是各向同性的，他们对载荷的变形反应取决于材料的分布方向，这种材料变形行为是材料内部超微小结构造成的结果，也就是跟纤维束有关，最常见的是胶原纤维、弹性纤维。弹性纤维主要存在于韧带和脉管壁，它与胶原纤维共同存在，赋予组织以弹性和抗张能力。弹性纤维如同橡皮带一样，它的长度能够伸展到正常长度的几倍，当收缩时又能恢复到原始长度。组织的弹性则是通过改变散布在弹性纤维中胶原的数量来控制。

图（4）肌肉结构

考察肌肉材料，他是各向异性材料的特例，他只有一个变形方向，我们把肌肉这类软体材料定义为横向各向同性材料，横向各向同性是指的是在垂直于轴线的任何方向上弹性相同，有5个独立的弹性常数。

肌肉具备超弹性，也叫Green弹性材料，是一种理想弹性材料的本构模型，其应力-应变关系源自应变能函数，需要注意的是超弹性材料是柯西弹性材料的特例。对于许多材料，线弹性模型（胡克方程，x = -kw）不能准确地描述观察到的材料行为。这类材料最常见的例子是橡胶，其应力-应变关系可定义为非线性弹性、各向同性、不可压缩且通常与应变率无关，超弹性为此类材料的应力-应变行为建模提供了理论框架。

图（5） Mooney 超弹性模型首次发表于 1940年

从20世纪40年代开始，弹性材料的有限变形理论得到了很大的发展。所得出的有意义理论结果和许多实验验证的结果被广泛地应用于描述橡胶类材料的物理行为和工程领域。弹性材料有限变形的数学理论的本质是非线性的，在理论和应用中所遇到的数学方面的困难是可想而知的。近年来，有许多力学工作者用数学和数值分析方法对许多有限弹性理论的非线性问题和技术难题进行了研究。在研究中，人们越来越认识到建立反映橡胶类材料变形特征本构模型是解决问题的核心，选取简单而适用的描述橡胶类材料的本构模型，决定着数值分析方法的精确性和得出结果的可用性。

1940年，在Isotropic、Isometric假设下，Melvin Mooney得到了适用于橡胶材料的应变能函数（如图5）（力学的论文都挺古老的）。1948年Rivlin从唯象学的角度出发对Mooney的工作做了改良，得出了Mooney-Rivlin模型，我们使用Mooney-Rivlin模型定义肌肉的应变能函数，但是因为肌肉的应变沿着纤维的方向伸缩，我们将其带入Mooney-Rivlin，加入一个矢量场用来表示应变方向。

标准的不可压缩Mooney-Rivlin材料的应变能函数如下：

公式1

$W = C_{1}(\bar{I_{1}}-3)+C_{2}(\bar{I_{2}}-3)$

其中 $C_{1}$ 及 $C_{2}$ 是材料参数，根据实验数据得来。 $\bar{I_{1}}$ 及 $\bar{I_{2}}$ 是公式 $\bar{B}=(det B)x^{\frac{-1}{b}}B$ 的一阶、二阶不变量。公式1之所以减去3的原因是，当无变形时，应变张量为单位张量 $I$ ，由公式1知其第一不变量和第二不变量为： $I_{1}=\delta ii = 3,I_{2}=\frac{1}{2}(\delta_{ii}^{2}-\delta_{ii}) = 3$ 。

公式2

$\bar{I_{1}}=J\frac{-2}{3}I_{1}$ , $I_{1}=\lambda_{2}^{1}+\lambda_{2}^{2}+\lambda_{2}^{3}$

$\bar{I_{1}}=J\frac{-4}{3}I_{2}$ , $I_{2}=\lambda_{2}^{1}\lambda_{2}^{2}+\lambda_{2}^{2}\lambda_{2}^{3}+\lambda_{2}^{3}\lambda_{2}^{1}$

其中应变梯度 F 及 $J = det(F)=\lambda_{1}\lambda_{2}\lambda_{3}$ ， $J=1$ 为不可压缩材料。

图（6）矢量场，粉色矩形是各向同性面

对于横向各向同性材料，一点的应变和应力状态除了和局部的变形梯度张量 $F$ 有关外，还和Figure 6 中所示的纤维方向有关。肌肉沿着纤维的方向变形，具备对称性，因此可以用一个矢量场 $\alpha^{0}(X)$ 去定义它，假设纤维伸长比为 $\lambda$ ，则其值可由应变梯度跟矢量场得到：

$\lambda\alpha=F\cdot\alpha^{0}$ ，由于 $\alpha^{0}$ 是一个单位矢量，因此：

公式3

$\lambda^{2}\alpha\cdot\alpha = \lambda^{2} = \alpha^{0} \cdot F^{T} \cdot F \cdot \alpha = \alpha^{0} \cdot C \cdot \alpha^{0}$

上式表明，纤维的伸长比不但依赖于参考构型中纤维的单位适量 $\alpha^{0}$ ，而且还和Right Cauchy-Green 变形张量 $C$ 有关，到此，横向各向同性的定义已完成（Transverse Isotropy）。

修改后的应变能函数 $\bar{W}$ 的一般形式如下：

公式4

$\bar{W}(X,C,\alpha)=\bar{W}(X,QCQ^{T},Q\alpha^{0}\otimes \alpha^{0} \otimes Q^{T} ) \forall(Q,C)\in Orth^{+} \times Sym^{+}$

横向各项同性超弹性体的行为由如下五个应变度量的不变量来确定：

公式5

$I_{1} = trC,I_{2}=\frac{1}{2}[(trC)^{2}-trC^{2}],I_{3}=detC = J^{2}$

$I_{4}=\alpha^{0} \cdot C \cdot \alpha^{0} , I_{5} = \alpha^{0} \cdot C^{2} \cdot \alpha^{0}$

其中 $I_{1},I_{2},I_{3}$ 是Right Cauchy Stretch Tensor，它表达了各向同性材料的不变量。

$I_{4}$ 是沿肌肉纤维方向两倍的拉伸强度， $I_{5}$ 定义了肌肉纤维束的各向异性特性，到此得出应变能函数 $\bar{W}$ 的最终形式：

公式6

$\bar{W}(X,C,\alpha^{0})=W(X,I_{1}(C),I_{2}(C),I_{3}(C),I_{4}(C,\alpha^{0}),I_{5}(C,\alpha^{0})) + \frac{1}{2}p(I_{3}-1)$

$p$ 是Lagrange算子，实际上是静水压强。

材料进阶（应力方程，弹性方程的推导）

图（7）分布在四面体各个面上的力

上一节得到了肌肉在受力后的变形规律，这一节要解决力的分布问题，在某个面上的应力是一个平均值，取决于微元的尺寸，尺寸越小精度越高。这里引入两个应力张量： $S = 2^{nd}$ Piola-Kirchhoff Stress Tensor， $\sigma$ = Cauchy Stress Tensor，前者表示变形前的应力情况，后者表示变形后的应力情况。它们之间的关系如下：

$\sigma^{ij} = \frac{1}{J}(FSF^{T})^{ij}$

$S^{ij} = J(F^{-1}\sigma^{-T})^{ij}$

不可压缩的超弹性材料的应力可以对其应变能泛函求偏导得到：

$S = 2\frac{\partial W}{\partial C} + pC^{-1}$

对公式6求偏导：

$S = 2\sum_{\alpha = 1}^{5}{\frac{\partial W}{\partial I_{\alpha}} \frac{\partial I_{\alpha}}{\partial C} } + pC^{-1}$

其中 $\frac{\partial I_{\alpha}}{\partial C}$ 根据公式5得出：

$\frac{\partial I_{1}}{\partial C} = 1 , \frac{\partial I_{2}}{\partial C} = I_{1} 1-C,\frac{\partial I_{3}}{\partial C} = I_{2}1-I_{1}C+C^{2}=I_{3}C^{-1},$

$\frac{\partial I_{4}}{\partial C} = \alpha^{0}\otimes \alpha^{0},\frac{\partial I_{5}}{\partial C}= \alpha^{0}\otimes C\cdot\alpha^{0} + \alpha^{0} \cdot C \otimes \alpha^{0}$

得到S = $2^{nd}$ Piola-Kirchhoff Stress Tensor， $\sigma$ = Cauchy Stress Tensor两个张量后，最后我们来定义弹性力学中的四阶张量，也就是本构模型里面的模量。这是因为对于一个质点，我们给定任意的二阶张量应变，能得到某种二阶张量应力，那么什么东西点乘二阶张量等于二阶张量呢？就是四阶张量了。详情可参考：[6]，四阶张量C：

$C=4 \{ (W_{11}+2W_{12}I_{1}+W_{2}+W_{22}I_{2}^{1})B\otimes B$

$-(W_{12}+W_{22}I_{1})(B \otimes B^{2} + B^{2} \otimes B) + W_{22}(B^{2} \otimes B^{2}) - W_{2}I_{B}$

$+(W_{14}+W_{24}I_{1}+W_{5})I_{4}(B \otimes \alpha \otimes \alpha + \alpha \otimes \alpha \otimes B)$

$+(W_{15}+W_{25}I_{1})(B \otimes [\varphi * \frac{\partial I_{5}}{\partial C} ] + [\varphi * \frac{\partial I_{5}}{\partial C} ] \otimes B)$

$-W_{24}I_{4}(B^{2}\otimes \alpha \otimes \alpha + \alpha \otimes \alpha \otimes B^{2})$

$-W_{25}(B^{2} \otimes [\varphi * \frac{\partial I_{5}}{\partial C} ] +[\varphi * \frac{\partial I_{5}}{\partial C} ] \otimes B^{2} )$

$-W_{44}I_{4}^{2}(\alpha \otimes\alpha \otimes\alpha \otimes\alpha) + W_{45}I_{4}(\alpha \otimes\alpha \otimes [\varphi * \frac{\partial I_{5}}{\partial C} ] + [\varphi * \frac{\partial I_{5}}{\partial C} ] \otimes \alpha \otimes \alpha)$

$+W_{55}( [\varphi * \frac{\partial I_{5}}{\partial C} ] \otimes [\varphi * \frac{\partial I_{5}}{\partial C} ]) \} + pI$

到此，不可压缩、超弹性、横向各向异性的应变能函数、应力张量、弹性张量的建模工作已经完成。接下来继续深入，为得出这些偏微分方程的数值解做准备。

有限单元法（Finite Element Method）

上一节定义的非线性偏微分方程（PDE）是非常复杂的，因为PDE有无穷多个解，所以在解决具体问题的时候，需要从中选取所需要的解，因此，还必须知道附加条件。

偏微分方程常用的解法有变分法（Caculus of Variations）和有限差分法（Finite Difference Method），变分法是把定解问题转化成变分问题，再求变分问题的近似解。有限元方法的基础是变分原理，我们用来模拟软组织，那其基本求解思想是把软组织（计算域Domain）划分为有限个互不重叠的单元（Element）如图9所示，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。

废话不多说，有限元的本质是什么？我们来看一个式子：

$\frac{d}{dx}(AE\frac{du}{dx})+b=0 ,on 0<x<1$ ，公式（7）

它有无穷多个解。

另外，它还出现了未知函数的导数。

所以求积分：

$\int_{0}^{1} \frac{d}{dx}(AE\frac{du}{dx})+b=0$ ，on $0<x<1$ 公式（8）

得先找到最优函数 $y(x)$ ，然后把它离散，最终就得到了数值解，就完成了全过程。
公式（7）是非常复杂的，我们要把它化简，过程是：

1.把二阶导数简化成一阶导数；

2.得到函数的变分形式，以便计算原函数的极值；

这个过程得到的就是原函数的 Weak Form，为什么叫弱形式呢？因为它的平滑度降低了。

$(wA\bar{t})_{x=0} + \int_{0}^{1}wbdx$ ，公式（9）

其中是 $\bar t$ nutrual boundary condition。

$wb$ 包含了essential boundary condition。

为什么加条件？方程里面全是未知数，选几个为数值预设几个值，解算就方便了，其实就是限定一下解的范围，这个在数学上，边界条件叫做定解条件。偏微分方程本身是表达同一类物理现象的共性，是作为解决问题的依据，定解条件却反映出具体问题的个性，它提出了问题的具体情况。方程和定解条件合而为一体，就叫做定解问题。

假设只有2个节点构成的单元，公式（9）离散之后：

$\sum_{e}^{N_{elment}} W^{eT} \{ { \int_{x_{1}^{e}}^{x_{2}^{e}} B^{eT} } A^{e}E^{e}B^{e}dxd^{e} - \int_{x_{1}^{e}}^{x_{2}^{e}} N^{eT}bdx -(N^{eT}A^{e}\bar t)_{x=0} \} = 0$

左边第一个积分是Stiffness Matrix $K^{e}$ ，第二个积分式是external Force $F_{\Omega^{e}}$ 。

他们的关系是：

$w^{T}(Kd-F) = 0$

最终肌肉顶点的 Displacements：

$d = W^{-1T}(W^{-1T}W^{-K})(W^{-1T}W^{F} )$

然后我们就可以解这个线性方程组了，有限元弱形式到离散化的推导过程就不给出了，以后有机会写文章介绍一下这个过程。

最后结语

图（11），吴老师一篇200多页的论文，已在Appiled Mechanics出版成书

工业软件，尤其是FEM，国内学术界有一个共识：不缺搞理论的，缺的是能把理论变成代码的，因此懂理论的不太擅长写代码，写代码的又不懂理论，就导致一些成果非常难落地。在此之前吴老师的一些研究都是找老外写代码，国内能写计算力学代码的人不多，无论学术界还是工业界大多数都用国外现成的软件做实验。另外，刷论文是唯一学术指标，体制内并不鼓励老师出来搞工程，甚至认为老师搞这些是不务正业。种种一些情况，导致国内FEM软件一直没有机会发展，我们做的是补齐这个短板，把一些领先的理论成果变成可用的软件产品，随着工程化的推进，技术慢慢形成了规范，就可以传承了，人才的培养跟得上，这块就可以良性发展。

References：

[1] A note on the elasticity of the soft biological tissue.J Biomechanics,5:309-311,1972.

[2] YC Fung.Elasticity of soft tissues in simple elongation.Am J Physiol,213:1532-1544,1967.

[3] RC Haut and RW Little.A constitutive equation for collagen fibers.J Biomechanics,5:421-430,1972.

[4] RB Jenkins and RW Little.A constitutive equation for parallel-fibered elastic tissue.J Biomechanics,7:397-402,1974.

[5] DR Veronda and R Westmann.Mechanical characterization of skin-finite doformations.J Biomechanics,3:111-124,1970.

[6] Mathematical Foundations of Elasticity (Dover Civil and Mechanical Engineering)

[7] 基本载荷作用下橡胶类材料的超弹性力学性能分析.