一、约束8-constraint.scn?xml version1.0? !-- Adding boundary conditions -- Node nameroot dt0.01 gravity0 0 0 Node nameplugins RequiredPlugin nameSofa.Component.AnimationLoop/ RequiredPlugin nameSofa.Component.Constraint.Projective/ RequiredPlugin nameSofa.Component.IO.Mesh/ RequiredPlugin nameSofa.Component.LinearSolver.Iterative/ RequiredPlugin nameSofa.Component.Mapping.NonLinear/ RequiredPlugin nameSofa.Component.Mass/ RequiredPlugin nameSofa.Component.MechanicalLoad/ RequiredPlugin nameSofa.Component.ODESolver.Backward/ RequiredPlugin nameSofa.Component.SolidMechanics.FEM.Elastic/ RequiredPlugin nameSofa.Component.StateContainer/ RequiredPlugin nameSofa.Component.Topology.Container.Dynamic/ RequiredPlugin nameSofa.Component.Visual/ RequiredPlugin nameSofa.GL.Component.Rendering3D/ /Node DefaultAnimationLoop computeBoundingBoxfalse/ VisualStyle displayFlagsshowForceFields/ !-- 粗体积网格用于力学计算 -- MeshGmshLoader namemeshLoaderCoarse filenamemesh/liver.msh/ !-- 细表面网格用于视觉渲染 -- MeshOBJLoader namemeshLoaderFine filenamemesh/liver-smooth.obj/ Node nameLiver EulerImplicitSolver/ CGLinearSolver iterations200 tolerance1e-09 threshold1e-09/ TetrahedronSetTopologyContainer nametopo src../meshLoaderCoarse/ TetrahedronSetGeometryAlgorithms templateVec3d nameGeomAlgo/ MechanicalObject templateVec3d nameMechanicalModel showObject1 showObjectScale3/ TetrahedronFEMForceField nameFEM youngModulus1000 poissonRatio0.4 methodlarge/ MeshMatrixMass massDensity1 topologytopo/ !-- 老师编辑器画面中增大到了100使变形更明显 -- ConstantForceField totalForce100 0 0/ !-- 节点1和50完全固定 -- FixedConstraint indices1 50/ !-- 节点3x方向自由y和z方向固定 -- PartialFixedConstraint indices3 fixedDirections0 1 1/ Node nameFineVisualModel OglModel nameVisualModel src../../meshLoaderFine/ BarycentricMapping nameVMapping input../MechanicalModel outputVisualModel/ /Node /Node /Node1.boundary conditions边界条件现在这个 SOFA 对象是 completely free to move也就是完全自由运动。只要有外力它整体就会被推走。但真实仿真里我们经常需要限制某些运动。例如某些节点固定在墙上、某些节点只能沿某个方向滑动、某个区域被夹具固定、某个边界不能动。这些就是boundary conditions边界条件。它们的作用是限制物体的自由运动。在 SOFA 中一种常见方式是用projective constraints投影约束其中一个最简单的就是FixedConstraint indices1 3 50 /SOFA 中还有其他类型的约束例如把一个点连接到另一个点、两个物体之间建立接触关系、两个自由度必须满足某种方程……这些约束不一定用FixedConstraint这种投影方式处理。有些约束需要另一套求解方法并且需要另一种 AnimationLoop后面还会再回来讲这部分。2.Projective Constraint 投影约束2.1.FixedConstraint 是一种投影约束FixedConstraint indices1 3 50 /这行代码 FixedConstraint indices1 3 50 / 意思是把编号为 1、3、50 的节点固定住。更具体一点这些节点的全部自由度都被固定。因为当前MechanicalObject是MechanicalObject templateVec3d ... /每个Vec3d节点有三个自由度 x,y,z完全固定表示这三个方向都不能运动​0,​0,02.2. 为什么叫“投影”约束按照一个时间步解释其工作过程。假设当前时间步是,大致过程是时间步开始 -- 把受约束节点的速度投影为0 -- 求解器计算新的位置和速度 -- 时间步结束 -- 再次把受约束节点投影回合法状态在时间步开始把节点1、3、50的速度投影为零完成求解、得到新位置和新速度后再次投影为零从而确保这些节点不会移动。3.检查被固定的节点运行场景后打开Behavior Models的显示可以看到受约束的三个节点被突出显示。这三个节点正是1、3、50。代码中的indices指向的是力学网格的节点编号而不是视觉网格上的三个任意点。TetrahedronSetGeometryAlgorithms templateVec3d nameGeomAlgo/用 GeometryAlgorithms 显示节点编号前面讲过它是一个几何工具包可以进行体积计算长度、面积等几何计算几何信息可视化。这里展示了它的另一个功能在界面中显示各个网格点的编号。4.固定节点后为什么物体仍然能变形现在作用在这个 liver 上的外力是 ConstantForceField!-- ConstantForceField totalForce1 0 0/ 为了让变形更明显又改成了 -- ConstantForceField totalForce100 0 0/这个 liver 已经具有MeshMatrixMass → 质量TetrahedronFEMForceField → 线弹性材料行为ConstantForceField → 外部节点力FixedConstraint → 限制部分节点运动ConstantForceField 是节点力( nodal force )固定三个点不等于整个物体变成刚体所以整个物体仍然可以产生弹性变形。5.PartialFixedConstraint 部分固定约束FixedConstraint会把一个节点的所有方向都锁死x 固定 y 固定 z 固定但有些物体需要像滑轨一样允许沿 x 方向移动禁止沿 y、z 方向移动这时使用PartialFixedConstraint/开始时三个节点全部固定FixedConstraint indices1 3 50/为了让节点3只受到部分限制需要先把节点3从完全固定中移除然后单独添加FixedConstraint indices1 50/ PartialFixedConstraint indices3 fixedDirections0 1 1/fixedDirections的三个数字依次对应x y z其中 0 不固定允许运动 1 固定禁止运动。这也是一种投影约束只不过它只投影部分自由度。二、重力9-gravity.scn?xml version1.0? !-- Gravity applied to a deformable object -- Node nameroot dt0.01 gravity0 -1 0 Node nameplugins RequiredPlugin nameSofa.Component.AnimationLoop/ RequiredPlugin nameSofa.Component.Constraint.Projective/ RequiredPlugin nameSofa.Component.IO.Mesh/ RequiredPlugin nameSofa.Component.LinearSolver.Iterative/ RequiredPlugin nameSofa.Component.Mapping.NonLinear/ RequiredPlugin nameSofa.Component.Mass/ RequiredPlugin nameSofa.Component.ODESolver.Backward/ RequiredPlugin nameSofa.Component.SolidMechanics.FEM.Elastic/ RequiredPlugin nameSofa.Component.StateContainer/ RequiredPlugin nameSofa.Component.Topology.Container.Dynamic/ RequiredPlugin nameSofa.Component.Visual/ RequiredPlugin nameSofa.GL.Component.Rendering3D/ /Node DefaultAnimationLoop computeBoundingBoxfalse/ VisualStyle displayFlagsshowForceFields/ MeshGmshLoader namemeshLoaderCoarse filenamemesh/liver.msh/ MeshOBJLoader namemeshLoaderFine filenamemesh/liver-smooth.obj/ Node nameLiver EulerImplicitSolver/ CGLinearSolver iterations200 tolerance1e-09 threshold1e-09/ TetrahedronSetTopologyContainer nametopo src../meshLoaderCoarse/ TetrahedronSetGeometryAlgorithms templateVec3d nameGeomAlgo/ MechanicalObject templateVec3d nameMechanicalModel showObject1 showObjectScale3/ TetrahedronFEMForceField nameFEM youngModulus1000 poissonRatio0.4 methodlarge/ !-- 重力必须通过质量转换成重力 -- MeshMatrixMass massDensity1 topologytopo/ FixedConstraint indices1 3 50/ Node nameFineVisualModel OglModel nameVisualModel src../../meshLoaderFine/ BarycentricMapping nameVMapping input../MechanicalModel outputVisualModel/ /Node /Node /Node1.Gravity接下来打开 9-gravity.scn场景根节点中有Node nameroot dt0.01 gravity0 -1 0三个数依次表示gx gy gzgravity0 -1 0 表示重力加速度沿 y 轴负方向。这里使用-1是为了演示。SOFA 不强制单位如果采用国际单位制地球重力通常会设置为接近: 0 -9.81 0 前提是网格、质量、时间和材料参数也全部使用一致单位gravity特点写的是加速度而不是直接写力必须通过质量得到重力质量越大对应的重力越大只有包含质量组件的对象才会受到重力根节点中的重力可以作用于其子树下所有有质量的力学对象。2.重力如何变成力重力加速度需要和质量结合才能得到重力从连续体角度更完整地说是把密度和重力加速度在物体体积内积分其中质量密度物体体积重力加速度​最终装配到节点上的重力。因此把重力称为volumetric external force field 体积外力场。它作用于整个有质量的物体而不是简单地手写一个节点力。三、Other physical model — FEM、弹簧模型、非均匀刚度与阻尼1.把 FEM 模型换成另一种物理模型TetrahedronFEMForceField nameFEM youngModulus1000 poissonRatio0.4/这个组件负责计算物体的线性弹性材料模型也就是材料本构关系。它把 liver 看成一个三维连续体通过四面体单元杨氏模量youngModulus泊松比poissonRatio计算物体受力后怎样变形、内部怎样产生恢复力。把这一行换成另一个组件MeshSpringForceField nameSprings tetrasStiffness10000/组件名字是 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并不一定只能是单个数值它可以是vectordouble 双精度数值向量情况一只有一个值youngModulus1000表示所有四面体都使用同一个杨氏模量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 个弹簧每一行对应一根具体弹簧。字段含义Index 1弹簧第一个端点的节点编号Index 2弹簧第二个端点的节点编号stiffnessKs弹簧刚度dampingKd弹簧阻尼L弹簧静止长度或参考长度刚度stiffness是材料力学中的名词定义为施力与所产生变形量的比值表示材料或结构抵抗变形的能力。公式记为其中表示刚度表示施力表示变形量。在国际单位制中刚度的单位为牛顿/米。一般应用于虎克定律作系统的振动分析。阻尼damping是指任何振动系统在振动中由于外界作用如流体阻力、摩擦力等和/或系统本身固有的原因引起的振动幅度逐渐下降的特性以及此一特性的量化表征。在实际振动中由于摩擦力总是存在的所以振动系统最初所获得的能量在振动过程中因阻力不断对系统做负功使得系统的能量不断减少振动的强度逐渐减弱振幅也就越来越小以至于最后的停止振动像这样的因系统的力学能由于摩擦及转化成内能逐渐减少振幅随时间而减弱振动称为阻尼振动。因此弹簧模型可以为不同弹簧设置不同KsFEM 模型可以为不同四面体设置不同youngModulus。两种模型都能做非均匀刚度只是控制对象不同3.3.弹簧组件里的阻尼MeshSpringForceField不仅能建立弹簧还能建立一个基础的 spring-damper model 弹簧—阻尼器模型。弹簧力主要取决于位移或伸长量阻尼力主要与相对速度有关。这里的边阻尼只是一个非常基础的阻尼模型不一定适合所有材料。“阻尼”可能有不同含义要区分两种“阻尼”物理阻尼如果你想模拟的是真实材料本身的黏性例如橡胶回弹有延迟生物组织有黏性软材料存在时间相关响应那么需要的是 Viscoelastic material 黏弹性材料这种阻尼属于真实材料本构模型的一部分。SOFA 中已经有人开发了用于定义黏弹性材料的插件。数值阻尼另一种情况是你并不一定要模拟真实黏性只是想让仿真更稳定更快收敛减少高频振荡避免求解结果一直抖动。这时使用的是 Numerical damping 数值阻尼比如 Rayleigh damping 瑞利阻尼数值阻尼的一种rayleighStiffness是 Rayleigh 数值阻尼用于帮助仿真稳定不一定代表真实材料的黏性。数值阻尼是为了数值稳定而添加的不一定对应真实材料中的黏性。甚至可以把它称为某种“fake damping”意思不是它没用而是它主要是数值计算手段不一定代表真实物理材料。