8-SOFA_约束、重力与材料模型补充

一、约束

8-constraint.scn

<?xml version="1.0"?> <!-- Adding boundary conditions --> <Node name="root" dt="0.01" gravity="0 0 0"> <Node name="plugins"> <RequiredPlugin name="Sofa.Component.AnimationLoop"/> <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <RequiredPlugin name="Sofa.Component.Mapping.NonLinear"/> <RequiredPlugin name="Sofa.Component.Mass"/> <RequiredPlugin name="Sofa.Component.MechanicalLoad"/> <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <RequiredPlugin name="Sofa.Component.StateContainer"/> <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <RequiredPlugin name="Sofa.Component.Visual"/> <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> </Node> <DefaultAnimationLoop computeBoundingBox="false"/> <VisualStyle displayFlags="showForceFields"/> <!-- 粗体积网格：用于力学计算 --> <MeshGmshLoader name="meshLoaderCoarse" filename="mesh/liver.msh"/> <!-- 细表面网格：用于视觉渲染 --> <MeshOBJLoader name="meshLoaderFine" filename="mesh/liver-smooth.obj"/> <Node name="Liver"> <EulerImplicitSolver/> <CGLinearSolver iterations="200" tolerance="1e-09" threshold="1e-09"/> <TetrahedronSetTopologyContainer name="topo" src="@../meshLoaderCoarse"/> <TetrahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo"/> <MechanicalObject template="Vec3d" name="MechanicalModel" showObject="1" showObjectScale="3"/> <TetrahedronFEMForceField name="FEM" youngModulus="1000" poissonRatio="0.4" method="large"/> <MeshMatrixMass massDensity="1" topology="@topo"/> <!-- 老师编辑器画面中增大到了100，使变形更明显 --> <ConstantForceField totalForce="100 0 0"/> <!-- 节点1和50完全固定 --> <FixedConstraint indices="1 50"/> <!-- 节点3：x方向自由，y和z方向固定 --> <PartialFixedConstraint indices="3" fixedDirections="0 1 1"/> <Node name="FineVisualModel"> <OglModel name="VisualModel" src="@../../meshLoaderFine"/> <BarycentricMapping name="VMapping" input="@../MechanicalModel" output="@VisualModel"/> </Node> </Node> </Node>

1.boundary conditions，边界条件

现在这个 SOFA 对象是 completely free to move，也就是完全自由运动。只要有外力，它整体就会被推走。

但真实仿真里，我们经常需要限制某些运动。例如：某些节点固定在墙上、某些节点只能沿某个方向滑动、某个区域被夹具固定、某个边界不能动。
这些就是boundary conditions，边界条件。它们的作用是限制物体的自由运动。

在 SOFA 中，一种常见方式是用：projective constraints，投影约束

其中一个最简单的就是：

<FixedConstraint indices="1 3 50" />

SOFA 中还有其他类型的约束，例如：把一个点连接到另一个点、两个物体之间建立接触关系、两个自由度必须满足某种方程……这些约束不一定用FixedConstraint这种投影方式处理。有些约束需要另一套求解方法，并且需要另一种 AnimationLoop，后面还会再回来讲这部分。

2.Projective Constraint 投影约束

2.1.FixedConstraint 是一种投影约束

<FixedConstraint indices="1 3 50" />

这行代码 <FixedConstraint indices="1 3 50" /> 意思是：把编号为 1、3、50 的节点固定住。

更具体一点：这些节点的全部自由度都被固定。因为当前MechanicalObject是：

<MechanicalObject template="Vec3d" ... />

每个Vec3d节点有三个自由度 x,y,z

完全固定表示这三个方向都不能运动：​=0,​=0,=0

2.2. 为什么叫“投影”约束

按照一个时间步解释其工作过程。

假设当前时间步是,

大致过程是：时间步开始 --> 把受约束节点的速度投影为0 --> 求解器计算新的位置和速度 --> 时间步结束 --> 再次把受约束节点投影回合法状态

在时间步开始，把节点1、3、50的速度投影为零；完成求解、得到新位置和新速度后，再次投影为零，从而确保这些节点不会移动。

3.检查被固定的节点

运行场景后，打开Behavior Models的显示，可以看到受约束的三个节点被突出显示。这三个节点正是：1、3、50。
代码中的indices指向的是力学网格的节点编号，而不是视觉网格上的三个任意点。

<TetrahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo"/>

用 GeometryAlgorithms 显示节点编号：
前面讲过，它是一个几何工具包，可以进行体积计算；长度、面积等几何计算；几何信息可视化。
这里展示了它的另一个功能：在界面中显示各个网格点的编号。

4.固定节点后为什么物体仍然能变形

现在作用在这个 liver 上的外力是 ConstantForceField

<!-- <ConstantForceField totalForce="1 0 0"/> 为了让变形更明显，又改成了 --> <ConstantForceField totalForce="100 0 0"/>

这个 liver 已经具有：
MeshMatrixMass → 质量
TetrahedronFEMForceField → 线弹性材料行为
ConstantForceField → 外部节点力
FixedConstraint → 限制部分节点运动

ConstantForceField 是节点力( nodal force )，固定三个点不等于整个物体变成刚体，所以整个物体仍然可以产生弹性变形。

5.PartialFixedConstraint 部分固定约束

FixedConstraint会把一个节点的所有方向都锁死：x 固定 y 固定 z 固定
但有些物体需要像滑轨一样：允许沿 x 方向移动；禁止沿 y、z 方向移动

这时使用<PartialFixedConstraint/>

开始时三个节点全部固定：

<FixedConstraint indices="1 3 50"/>

为了让节点3只受到部分限制，需要先把节点3从完全固定中移除，然后单独添加：

<FixedConstraint indices="1 50"/> <PartialFixedConstraint indices="3" fixedDirections="0 1 1"/>

fixedDirections的三个数字依次对应x y z，其中 0 = 不固定，允许运动 ；1 = 固定，禁止运动。

这也是一种投影约束，只不过它只投影部分自由度。

二、重力

9-gravity.scn

<?xml version="1.0"?> <!-- Gravity applied to a deformable object --> <Node name="root" dt="0.01" gravity="0 -1 0"> <Node name="plugins"> <RequiredPlugin name="Sofa.Component.AnimationLoop"/> <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <RequiredPlugin name="Sofa.Component.IO.Mesh"/> <RequiredPlugin name="Sofa.Component.LinearSolver.Iterative"/> <RequiredPlugin name="Sofa.Component.Mapping.NonLinear"/> <RequiredPlugin name="Sofa.Component.Mass"/> <RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <RequiredPlugin name="Sofa.Component.StateContainer"/> <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <RequiredPlugin name="Sofa.Component.Visual"/> <RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> </Node> <DefaultAnimationLoop computeBoundingBox="false"/> <VisualStyle displayFlags="showForceFields"/> <MeshGmshLoader name="meshLoaderCoarse" filename="mesh/liver.msh"/> <MeshOBJLoader name="meshLoaderFine" filename="mesh/liver-smooth.obj"/> <Node name="Liver"> <EulerImplicitSolver/> <CGLinearSolver iterations="200" tolerance="1e-09" threshold="1e-09"/> <TetrahedronSetTopologyContainer name="topo" src="@../meshLoaderCoarse"/> <TetrahedronSetGeometryAlgorithms template="Vec3d" name="GeomAlgo"/> <MechanicalObject template="Vec3d" name="MechanicalModel" showObject="1" showObjectScale="3"/> <TetrahedronFEMForceField name="FEM" youngModulus="1000" poissonRatio="0.4" method="large"/> <!-- 重力必须通过质量转换成重力 --> <MeshMatrixMass massDensity="1" topology="@topo"/> <FixedConstraint indices="1 3 50"/> <Node name="FineVisualModel"> <OglModel name="VisualModel" src="@../../meshLoaderFine"/> <BarycentricMapping name="VMapping" input="@../MechanicalModel" output="@VisualModel"/> </Node> </Node> </Node>

1.Gravity

接下来打开 9-gravity.scn

场景根节点中有：

<Node name="root" dt="0.01" gravity="0 -1 0">

三个数依次表示：gx gy gz
gravity="0 -1 0" 表示重力加速度沿 y 轴负方向。
这里使用-1是为了演示。SOFA 不强制单位，如果采用国际单位制，地球重力通常会设置为接近: 0 -9.81 0 （前提是网格、质量、时间和材料参数也全部使用一致单位）

gravity特点：

• 写的是加速度，而不是直接写力；
• 必须通过质量得到重力；
• 质量越大，对应的重力越大；
• 只有包含质量组件的对象才会受到重力；
• 根节点中的重力可以作用于其子树下所有有质量的力学对象。

2.重力如何变成力

重力加速度需要和质量结合，才能得到重力：

从连续体角度，更完整地说是把密度和重力加速度在物体体积内积分：

其中：

• ：质量密度；
• ：物体体积；
• ：重力加速度；
• ​：最终装配到节点上的重力。

因此把重力称为：volumetric external force field 体积外力场。它作用于整个有质量的物体，而不是简单地手写一个节点力。

三、Other physical model — FEM、弹簧模型、非均匀刚度与阻尼

1.把 FEM 模型换成另一种物理模型

<TetrahedronFEMForceField name="FEM" youngModulus="1000" poissonRatio="0.4"/>

这个组件负责计算物体的线性弹性材料模型，也就是材料本构关系。

它把 liver 看成一个三维连续体，通过：

• 四面体单元；
• 杨氏模量youngModulus；
• 泊松比poissonRatio；

计算物体受力后怎样变形、内部怎样产生恢复力。

把这一行换成另一个组件：

<MeshSpringForceField name="Springs" tetrasStiffness="10000"/>

组件名字是 Mesh Spring Force Field 网格弹簧力场，它是一个弹簧模型。

原来： TetrahedronFEMForceField ——> 把物体当成连续弹性材料

现在：MeshSpringForceField ——> 把网格节点之间的边看成弹簧

虽然场景中的 Loader、Topology、MechanicalObject、Mass、FixedConstraint、VisualModel、Mapping 都没有改变，但只替换这一行，物体的变形计算方法已经完全改变。

2.弹簧模型

FEM 模型认为这是一个连续物体，需要在空间中、在每个单元内部进行积分，并使用杨氏模量和泊松比等材料参数。

弹簧模型则更简单：不再仔细描述连续材料，直接在网格的边上放弹簧。

例如网格中有三个节点：

A ●────────● B \ / \ / \ / ● C

边：

A—B A—C B—C

都可以被看成弹簧。

因此：
节点 = 弹簧连接的位置
边 = 弹簧本身

这是一个非常基础的线性模型，不是像 FEM 那样对三维连续材料进行体积积分。

3.设置不同区域刚度

3.1.FEM 设置不同区域刚度

查看TetrahedronFEMForceField的属性。youngModulus并不一定只能是单个数值，它可以是：
vector<double> 双精度数值向量

情况一：只有一个值

youngModulus="1000"

表示所有四面体都使用同一个杨氏模量：

tetra 0 → 1000 tetra 1 → 1000 tetra 2 → 1000 ……

这就是均匀材料。

情况二：提供多个值

也可以为不同四面体提供不同的杨氏模量：

tetra 0 → 500 tetra 1 → 500 tetra 2 → 2000 tetra 3 → 2000

这样同一个物体中就可以出现一部分软，一部分硬。

3.2.弹簧模型如何设置不同刚度

打开MeshSpringForceField的属性窗口。

截图中可以看到一个弹簧列表，有914 个弹簧，每一行对应一根具体弹簧。

刚度（stiffness）是材料力学中的名词，定义为施力与所产生变形量的比值，表示材料或结构抵抗变形的能力。公式记为

其中表示刚度，表示施力，表示变形量。在国际单位制中，刚度的单位为牛顿/米。一般应用于虎克定律作系统的振动分析。

阻尼（damping）是指任何振动系统在振动中，由于外界作用（如流体阻力、摩擦力等）和/或系统本身固有的原因引起的振动幅度逐渐下降的特性，以及此一特性的量化表征。

在实际振动中，由于摩擦力总是存在的，所以振动系统最初所获得的能量，在振动过程中因阻力不断对系统做负功，使得系统的能量不断减少，振动的强度逐渐减弱，振幅也就越来越小，以至于最后的停止振动，像这样的因系统的力学能，由于摩擦及转化成内能逐渐减少，振幅随时间而减弱振动，称为阻尼振动。

因此：

• 弹簧模型可以为不同弹簧设置不同Ks；
• FEM 模型可以为不同四面体设置不同youngModulus。

两种模型都能做非均匀刚度，只是控制对象不同

3.3.弹簧组件里的阻尼

MeshSpringForceField不仅能建立弹簧，还能建立一个基础的 spring-damper model 弹簧—阻尼器模型。
弹簧力主要取决于位移或伸长量，阻尼力主要与相对速度有关。

这里的边阻尼只是一个非常基础的阻尼模型，不一定适合所有材料。

“阻尼”可能有不同含义，要区分两种“阻尼”：

物理阻尼

如果你想模拟的是真实材料本身的黏性，例如：

• 橡胶回弹有延迟；
• 生物组织有黏性；
• 软材料存在时间相关响应；

那么需要的是 Viscoelastic material 黏弹性材料

这种阻尼属于真实材料本构模型的一部分。SOFA 中已经有人开发了用于定义黏弹性材料的插件。

数值阻尼

另一种情况是：你并不一定要模拟真实黏性，只是想让仿真：

• 更稳定；
• 更快收敛；
• 减少高频振荡；
• 避免求解结果一直抖动。

这时使用的是 Numerical damping 数值阻尼，比如 Rayleigh damping 瑞利阻尼（数值阻尼的一种）
rayleighStiffness是 Rayleigh 数值阻尼，用于帮助仿真稳定，不一定代表真实材料的黏性。

数值阻尼是为了数值稳定而添加的，不一定对应真实材料中的黏性。甚至可以把它称为某种“fake damping”，意思不是它没用，而是它主要是数值计算手段，不一定代表真实物理材料。