草庐IT

第三讲:刚体动力学模拟

离原春草 2023-04-11 原文

上一讲中介绍了单个block或者两个block的弹簧系统在受力情况下的运动过程,其中介绍了三种数值积分的方法:

  • 显式欧拉
  • 隐式欧拉
  • 辛欧拉

不过前面的内容中物理系统的自由度较少,在实际情况中我们遇到的情况比这个复杂的多,因此今天这里来介绍刚体在受力情况下的运动过程,刚体指的是在任何受力情况下都不会发生形变的物体,这是一种理想情况,实际情况中并不存在刚体这种物体。

不过这里值得一提的是,将物体看成刚体的假设实际上也会引入新的问题,我们后面会看到,刚体的计算复杂度比软体更高,这也是为什么现有的刚体模拟中,没有哪种方法是真正完全准确的原因。另外,作为一般性,我们会发现,有时候我们在解决问题中的假设,试图将问题进行简化,但实际上并没有得到很好的效果,这里的刚体假设是一个案例,另一个案例就是后面介绍的流体动力学,在流体中,我们给出的一个假设是流体不可压缩,而这个假设同时也加剧了流体计算的复杂度。

用数学公式来表示刚体,可以给出如下:

即对物体上的任意两点,其欧式距离都是一个固定的常量。

在上面的情况下,我们可以用一个很长的向量来表示刚体:

要计算这个很长的向量就会变得复杂,因此一个很直观的方法就是采用reduced coordinate(广义坐标)来表达,简单描述就是,由于任意两点之间的固定位置关系,因此可以用一些较少的自由度来完成对上述长向量的表达,后面我们会介绍,对于一个刚体来说,我们通常只需要六个自由度就能够完成表达。

而广义坐标方法的使用策略给出如下:

  1. 定义广义坐标
  2. ”immersion“映射,通过表达式将定义的广义坐标映射到full coordinate
  3. 根据表达式,就可以定义刚体的动能
  4. 根据上面的内容,就可以写出外力做功的表达式
  5. 最终输出一个与牛顿第二定律等价的广义坐标方程

根据上面的广义坐标方程,就可以利用上一讲的内容利用数值积分完成求解,上面策略中前两步称为kinematics,后面三步称为dynamics。

从上面的策略,我们看到,这里是从能量角度而非力的角度(因为在实际的情况中,力的情况已经十分复杂且不清楚了)出发推导出一个广义坐标方程,而这就是经典物理中的一个重要里程碑拉格朗日力学的最重要的思想——第一次从能量的角度来看待物理世界。这样的角度的好处是,能量是跟坐标系无关的,所以推导跟建模的过程就无需考虑广义坐标的具体形式是什么样子。

下面会介绍两块内容:

  1. kinematics在广义坐标系下是什么样的
  2. 时间不够,这里会跳过拉格朗日力学的具体推导过程,即一个经典的变分思想,直接给出牛顿第二定律在广义坐标下的形式

这一讲的重点在于外力,而游戏中的外力99%的情况下来自于外界的碰撞,因此这里会主要介绍如何解决碰撞。

1. 广义坐标使用策略细节

首先来定义一下计算的坐标系:

上图给出了一个简单的坐标系示意图,对于刚体上的任意一点,我们可以用一个3D向量来表示,而在这个表示方法下,我们可以定义一个刚体的质心(center of mass)为:

相当于对每一点的坐标进行了加权平均,而权重是当点质量的百分比。

这里来定义广义坐标系,而广义坐标系的原点,我们就选取在物体初始状态时的质心位置,即,这是大家常用的做法,在后面的计算中按照这种做法会使得计算得到一定的简化。

刚体在初始状态受力之后,位置会发生变化,我们定义初始状态的刚体上某点位置叫,而变化后的位置可以用如下公式计算得到:

其中指的是围绕质心旋转某个角度的旋转矩阵,也就是说,整个刚体的运动可以拆解成平移与旋转,平移通过质心的坐标变化来表示,而旋转都是绕着质心完成的。

刚体上每个质点的动能为,其中为质点的速度,而这个速度可以通过如下公式计算得到:

可以看到,质点的速度可以由质心的速度(上面等式右边的求和中的右边一项)与质点的旋转速度表示。

在2D的情况下,旋转矩阵是一个2x2的矩阵:

导数就可给出为:

3D坐标系下情况要复杂一点,是一个三维旋转矩阵,属于Lie Group(李群),Lie Group表示的某一类矩阵,而Lie algebra(李代数)是另外一个代数结构,表示另外一类矩阵,这里用A表示。

我们知道,三维旋转可以由旋转轴与旋转角度表示,旋转轴可以表示为:
,旋转角度这里就还是用表示。根据旋转轴跟旋转角度,我们可以构建出一个Lie Algebra的矩阵:

这是一个斜对称矩阵(Skew-symmetric matrix),即其转置等于矩阵本身取负:

注意,这里的A不是旋转矩阵,但是跟我们的旋转矩阵有一个很好的映射关系:,这里的矩阵指数计算跟标量指数的计算是一样的,不同的是矩阵指数是将指数应用到矩阵上的每个元素上。

假设某个质点在当前的旋转状态为,其旋转角速度为,角速度可以用两个量进行表达,即旋转轴与角速度大小,令,这个三维坐标可以看成是旋转坐标乘上旋转的角速度大小,那么经过一个短暂时间的作用之后,其旋转的角度就变成了,其对应的Lie algebra就等于:

那么在状态下,经过了时间dt之后的旋转状态为,那么R_t的导数就可以表示为:

进行泰勒展开,即,省去后面的子项,上面公式可以变成:

将A(dt)代入,就得到了:

而:

另外,需要说一下的是,

有了速度,我们就可以计算出动能:



又因为为初始质心坐标,即坐标原点,因此这个求和项结果为0,所以:

分别称为Translational Energy与Rotational Energy。用矩阵形式来表达:

这里是一个对角矩阵,矩阵中的每个元素都等于

这里的被称为3D moment of inertia(通常使用inertia的首字母表示,inertia是惯性,也可以看成是另一种形式的质量,这里的inertia是旋转时的质量),这里会有一个简单的推导过程,时间关系,这里就直接省掉了。

回到微分方程:

令动量,两边乘以m,我们有:

这是质点移动的微分方程,总共包含了6个自由度(各三个),如果我们要为刚体建立一个类似的表达式,我们就需要将前面表示动能的所有自变量总共12个自由度都考虑进去:

  • 质点位置 (后面考虑省略上面的箭头以简化书写)三个自由度
  • 旋转轴 也是三个自由度(虽然多了个 ,但是由于旋转轴方向长度为1,所以又减去了一个自由度,刚好还是3个),
  • 刚体线速度 三个自由度
  • 刚体角速度 三个自由度)

同样,这里用动量来取代速度:为线动量,为角动量,因此我们就可以得到:

这里来逐项解释一下:

  • 第一项表示质心的位移的导数,可以理解为刚体的线速度,因此可以转化为动量除以质量,这里的质量是整个刚体的质量
  • 第二项是角度的导数,对应的就是角速度,而角速度前面已经求解过,可以直接写入
  • 第三项动量的导数,也就是的导数,就是也就是质心的受力
  • 第四项角动量的导数,那就是质心的力矩了,力矩可以在网上搜到求解方法,这里就不做赘述。

这就是刚体下的牛顿第二定律。

2. 碰撞

游戏中最常见的物理模拟就是刚体的碰撞,碰撞主要有两个计算阶段:

  • 碰撞检测,判断是否发生碰撞,以及碰撞发生的位置
  • 受力计算,碰撞发生后会对产生碰撞的两个物体施加怎样的力和力矩

碰撞检测的方法在网上可以搜到大量资料,这里因为时间关系,这里就不做展开了,接下来的内容主要介绍如何计算受力,而这个问题到现在为止还没有一个完美的解决方案,现有的物理引擎提供了一系列的solver,允许用户进行参数调节,但是这些solver都有各自的缺陷。

两种求解方法:

  • Impulse-based方法:大概思路是通过冲量来完成对碰撞影响的模拟,简单来说,设碰撞前某点的速度为,碰撞后的速度为,在一个time step(通常是一个物理帧间隔)之内,碰撞对二者产生了一个冲量,经过这个冲量导致两者的速度发生了变化
  • Penalty-based方法,这是游戏引擎中常用的方法,其大致思路是在碰撞发生的时候为双方施加一个penalty,可以理解为在碰撞的时候,两者发生挤压,从而产生弹力,而弹力就会对应于弹性势能这里的h理解为两者碰撞的时候挤压的幅度(想象弹簧的形变幅度),而通过对势能求导我们可以得到对应的受力:

先来看下Penalty-based方法,在Penalty-based方法中:

  • 假设我们就用弹性势能作为碰撞时候产生的势能,那么我们就有:,通过这个受力来计算碰撞发生后两者的速度与位置,而这种方法的一个问题在于,如果发生碰撞时的速度过快,在给定的k值下,可能没有办法将两者弹开,而是可能会导致一个物体从另一个物体中穿过去的异常现象;
  • 如上图所示,前面的弹性势能对应的曲线用红线绘制,可以看到势能增长的幅度与挤压的幅度的二次方程正比,而如果我们能够给出一种penalty function(称之为barrier function)使得势能在接近某个阈值的时候迅速逼近无穷大就能够有效抑制前面的问题。不过这种算法理想是美好的,现实却很残酷,比如在数值计算中没有办法处理这种无穷大的逻辑

在Imulse-based方法中,有其中是碰撞点的法线(方向朝外),是当前的冲量。为了避免出现前面Penalty-based方法中物体穿过碰撞体的情况,这里通常会给碰撞施加一个约束,比如发生碰撞的时候的一个质点的位置为q,那么我们就有:

这里g(q)可以表示为当前点高于碰撞面往外方向的长度,举个最简单的例子,某个刚体与平面发生碰撞,我们可以有:

上面公式中表示的是碰撞面上的点,则是碰撞面上的法线。在Impulse-based方法中可以通过冲量很容易就能保证这种约束条件的达成。

那么如何求取这里的冲量呢(或者换句话说,如何计算输出的速度)?这里有很多模型来解决这个问题,最经典的做法是根据能量守恒定律,物体碰撞前跟碰撞后的能量保持不变:那么不考虑碰撞时的能量损耗,可以直接将输出速度的幅度设置为与输入速度一样,方向则可以根据反射的方式求得(当然,这是在碰撞没有考虑力矩作用的前提下的结论);实际情况下,我们为了保证模拟结果的精度,通常会考虑碰撞过程中的能量损耗(如摩擦导致的热能损耗等),这种时候我们会添加一个叫做restitution coefficient的系数来控制能量保留比例,如果为0,则表示能量完全损耗。

上面介绍的是单点碰撞的情况,但我们有时候需要处理多点碰撞(multiple contact)的问题,这种问题就要复杂得多,这里就先不讲了,看看后面课程会不会有相关的内容。

有关第三讲:刚体动力学模拟的更多相关文章

  1. ruby - 如何模拟 Net::HTTP::Post? - 2

    是的,我知道最好使用webmock,但我想知道如何在RSpec中模拟此方法:defmethod_to_testurl=URI.parseurireq=Net::HTTP::Post.newurl.pathres=Net::HTTP.start(url.host,url.port)do|http|http.requestreq,foo:1endresend这是RSpec:let(:uri){'http://example.com'}specify'HTTPcall'dohttp=mock:httpNet::HTTP.stub!(:start).and_yieldhttphttp.shou

  2. ruby-on-rails - 在这种情况下我如何模拟一个对象?没有明显的方法可以用模拟替换对象 - 2

    假设我在Store的模型中有这个非常简单的方法:defgeocode_addressloc=Store.geocode(address)self.lat=loc.latself.lng=loc.lngend如果我想编写一些不受地理编码服务影响的测试脚本,这些脚本可能已关闭、有限制或取决于我的互联网连接,我该如何模拟地理编码服务?如果我可以将地理编码对象传递到该方法中,那将很容易,但我不知道在这种情况下该怎么做。谢谢!特里斯坦 最佳答案 使用内置模拟和stub的rspecs,你可以做这样的事情:setupdo@subject=MyCl

  3. ruby - "public/protected/private"方法是如何实现的,我该如何模拟它? - 2

    在ruby中,你可以这样做:classThingpublicdeff1puts"f1"endprivatedeff2puts"f2"endpublicdeff3puts"f3"endprivatedeff4puts"f4"endend现在f1和f3是公共(public)的,f2和f4是私有(private)的。内部发生了什么,允许您调用一个类方法,然后更改方法定义?我怎样才能实现相同的功能(表面上是创建我自己的java之类的注释)例如...classThingfundeff1puts"hey"endnotfundeff2puts"hey"endendfun和notfun将更改以下函数定

  4. ruby - 在 RSpec 中 stub /模拟全局常量 - 2

    我有一个gem,它有一个根据Rails.env的不同行为的方法:defself.envifdefined?(Rails)Rails.envelsif...现在我想编写一个规范来测试这个代码路径。目前我是这样做的:Kernel.const_set(:Rails,nil)Rails.should_receive(:env).and_return('production')...没关系,只是感觉很丑。另一种方法是在spec_helper中声明:moduleRails;end而且效果也很好。但也许有更好的方法?理想情况下,这应该有效:rails=double('Rails')rails.sho

  5. ruby-on-rails - rspec 模拟对象属性赋值 - 2

    我有一个rspec模拟对象,一个值赋给了属性。我正在努力在我的rspec测试中满足这种期望。只是想知道语法是什么?代码:defcreate@new_campaign=AdCampaign.new(params[:new_campaign])@new_campaign.creationDate="#{Time.now.year}/#{Time.now.mon}/#{Time.now.day}"if@new_campaign.saveflash[:status]="Success"elseflash[:status]="Failed"endend测试it"shouldabletocreat

  6. ruby - 如何使用 rspec stub /模拟对命令行的调用? - 2

    我正在尝试测试命令行工具的输出。如何使用rspec来“伪造”命令行调用?执行以下操作不起作用:it"shouldcallthecommandlineandreturn'text'"do@p=Pig.new@p.should_receive(:run).with('my_command_line_tool_call').and_return('resulttext')end如何创建stub? 最佳答案 使用newmessageexpectationsyntax:规范/虚拟规范.rbrequire"dummy"describeDummy

  7. 建模分析 | 平面2R机器人(二连杆)运动学与动力学建模(附Matlab仿真) - 2

    目录0专栏介绍1平面2R机器人概述2运动学建模2.1正运动学模型2.2逆运动学模型2.3机器人运动学仿真3动力学建模3.1计算动能3.2势能计算与动力学方程3.3动力学仿真0专栏介绍?附C++/Python/Matlab全套代码?课程设计、毕业设计、创新竞赛必备!详细介绍全局规划(图搜索、采样法、智能算法等);局部规划(DWA、APF等);曲线优化(贝塞尔曲线、B样条曲线等)。?详情:图解自动驾驶中的运动规划(MotionPlanning),附几十种规划算法1平面2R机器人概述如图1所示为本文的研究本体——平面2R机器人。对参数进行如下定义:机器人广义坐标

  8. ruby - 接收 block 作为参数的模拟方法 - 2

    我有一个或多或少这样的场景classAdefinitialize(&block)b=B.new(&block)endend我正在对A类进行单元测试,我想知道B#new是否正在接收传递给A#new的block。我使用Mocha作为模拟框架。这可能吗? 最佳答案 我用Mocha和RSpec都试过了,虽然我可以通过测试,但行为不正确。从我的实验中,我得出结论,验证block是否已通过是不可能的。问题:为什么要传递一个block作为参数?block将用于什么目的?什么时候调用?也许这确实是您应该用类似的东西测试的行为:classBlockP

  9. ruby - 如何模拟 Fixnum 变量的整数溢出? - 2

    我目前正在将一种算法从Java转换为Ruby,但由于Ruby中缺少整数溢出,我遇到了一些障碍。假设我的值为2663860877,它大于最大整数2147483648。在Java中,它环绕,我应该得到-1631106419。我找到了这段代码,但它似乎不起作用:defforce_overflow(i)ifi2147483647i&0xffffffffelseiendend并且'ing变量不会像您期望的那样强制它为负。 最佳答案 假设32位整数具有二进制补码负数,这应该可行:defforce_overflow_signed(i)force_

  10. ruby - 单元测试 ruby​​ 命令行应用程序的代码 - 如何模拟/传递 ARGV - 2

    我有一个命令行应用程序,它使用thor来处理选项的解析。我想使用test-unit和/或minitest针对代码对命令行功能进行单元测试。我似乎无法弄清楚如何确保ARGV数组(通常会保存命令行中的选项)保存我的测试选项,以便它们可以根据代码进行测试。具体应用代码:#myapp/commands/build.rbrequire'thor'moduleMyappmoduleCommands#DefinebuildcommandsforMyAppcommandlineclassBuild:test_unit#Definesourcerootofapplicationdefself.sourc

随机推荐