3D Game Programming Design(六)：物理系统与碰撞

老实说物理方面的研究对我来说实际上是个很严峻的课题，简单的物理模拟尚可以用我仅剩的高中牛顿物理知识计算解决，但稍微复杂一点的物理模拟我就只能举手投降了。

我们先来看看这次的我们要做的小游戏

打靶游戏：

游戏内容要求：
1. 靶对象为 5 环，按环计分；
2. 箭对象，射中后要插在靶上
3. 游戏仅一轮，无限 trials
4. 添加一个风向和强度标志，提高难度
5. 射中后，箭体产生颤抖效果

前面三项要求颇为简单，靶分数用一个到靶心的radius的计算分数即可。插在靶上用碰撞体解决。无限trial等于没说。箭的下落可以用quaternion.fromtorotation。如果箭模型本身是分离的话，可能还需要使用Unity提供的joint脚本（实际上我很恐惧使用这个东西，因为我曾经用过它来控制我制作的一个链球，却始终无法让链球的部分在地面拖动时可以保持一个正常的物理表现）

问题在模拟真实的风向和箭体的颤抖效果，如果我们需要使用joint的话还需要了解如何控制类似弹簧和刚体弹簧的表现。刚好前段时间在图书馆里借了一本不少人推荐的游戏用物理书，干脆就利用这个机会一边看书一边实践吧。（可怜实时渲染第三版又被我插队，不过有句讲句，那玩意儿真的好难懂……刷完一整本shader入门精要只能算勉强入了门，自己来打代码需要不断地翻笔记）

五天后：

好吧，我放弃用数学方法在unity里模拟物理了，实在是太太太太太糟心了，还是直接用已经给了我的物理引擎去模拟比较靠谱……

简单做了一下，重心效果和抖动效果如图，其实和我之前看的书没啥大关系……分割线下的东西可以无视

Unity我找不到可以给模型刷密度的功能，所以干脆做了个小球当做箭头的重心，用fixed joint（固定铰链）和箭的模型连在一起，然后修改箭的阻力和球的阻力，好处是不需要自己处理重力和空气阻力，坏处是交给物理引擎模拟，连我都不知道这玩意儿的落点究竟在哪里……

抖动因为搞不懂物理上的原因，所以用感觉像是不倒翁一样的代码控制。每次update向随机方向转动，并且减少转动幅度直到time结束

shake.cs：

using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class shake : MonoBehaviour {

    public Transform transform;
    public float time = 0; //时长
    public float flood = 0; //幅度
    public float speed = 0; //速度

    private Vector2 dstpos;
    private Vector2 currentpos;
    private float curflood;
    private float curtime;

	void Start () {
        if (transform == null) {
            transform = this.GetComponent<Transform>();
        }

        currentpos = Vector2.zero;
        dstpos = currentpos;
        curflood = flood;
        curtime = time;
	}
	
	void Update () {
        if (curtime <= 0) return;
        curtime -= Time.deltaTime;

        float curspeed = speed * Time.deltaTime;
        curflood -= (Time.deltaTime / time) * flood; //减少振幅

        if (dstpos == currentpos) dstpos = Random.onUnitSphere * curflood; //如果停止，重新分配随机dst

        if (Vector2.Distance(dstpos, currentpos) > curspeed) 
            currentpos += (curspeed / Vector2.Distance(dstpos, currentpos)) * (dstpos - currentpos) ;
        else currentpos = dstpos;

        transform.rotation = Quaternion.Euler(new Vector3(currentpos.x, currentpos.y, 0) + transform.eulerAngles); //转动角度
	}
}

横向风力影响很简单，用addForceAtPosition即可

导致弓转动的风力影响我看到后面书有写，不过感觉挺难懂的，这个效果也可以了，先鸽着吧（我知道这一鸽很有可能就是永远）

————————————————————————分割线————————————————————————

几个比较重要但不熟悉的术语

重心（质心），大家很熟悉这个概念意味着什么，但只有很少人可以说出质心的特性——当任何外力施加在该点上时，都不会引起物体的转动。重心的计算方法是每个质块的相对于右手参考坐标系的一阶矩的和除以物体总质量，得到的是一个维度上的质心的坐标，你需要每个维度都做这样的计算，一阶矩的计算方法如下：先得到沿着某个特定坐标轴从原点到质块中心的距离，该距离与该质块的质量的乘积即为一阶矩
转动惯量，衡量物体对转动的阻力。计算方法是每个质块对每条坐标轴的二阶矩乘以转动惯量相对的坐标到质块中心的垂直距离的平方乘以质量。看不懂？没关系，因为这玩意儿太复杂，我们一般不会去计算他。我们有一条实用观测公式：给定物体对于穿过物体质心的中性轴的转动惯量为I0，那么对于平行于这条中性轴的一定距离的轴的转动惯量I计算公式为—— I = I0 + m*d^2，m是物体的质量，d是两条平行轴的垂直距离。我们通常会用近似简单物体代替非规则物体来进行近似求解，进一步，将复杂物体分成小组件，并利用某些组件的I0相对于其md^2对于整个物体的转动惯量可以忽略的简化计算。至于规则物体的转动惯性，你可以在动力学教材里找到数据。

下面我们举一个例子来示范下这两个参数是如何计算的。假设我们现在要计算一辆卡车，如图

三个物件的参数如下图

不难计算出总重量Wtotal = 19343N，默认9.8加速度处理，总质量Mtotal = 1972 kg

整体的重心位置是不难计算的，我们需要分别对X轴和Y轴计算一次一阶矩

Xcgbody = (17500 * 30.5 + 850 * 31.5 + 993 * 28 ) / 19343 = 30.42 m

Ycgbody = (17500 * 30.5 + 850 * 31 + 993 * 30.5) / 19343 = 30.52 m

重心坐标即为(30.42, 30.52)

我们下面计算转动惯量，因为我们不希望车辆侧翻，所以只计算三维Y轴的转动惯量。

第一步：计算每个组件相对于自己中性轴的局部转动惯量。我们把三个组件都近似看作方柱体，已知方柱体的转动惯量公式如下图

Icar = (17500 / 9.8 / 12) * ((1.8)^2 + (4.7)^2) = 3765.5

Idriver = (850 / 9.8 / 12) * ((0.9)^2 + (0.5)^2) = 7.7

Ifuel = (993 / 9.81 / 12) * ((0.5)^2 + (0.8)^2) = 8.9

第二步，我们需要计算到物体重心的中性轴的转动惯量，我们先计算组件重心到物件重心的距离的平方（过程省略）

dcar^2 = 0.01

ddriver^2 = 1.68

dfuel^2 = 5.86

第三步，代入计算平行轴定理 I = I0 + m*d^2（过程省略）

Icar = 3783.34

Idirver = 153.27(7.7的中性轴转动惯量占比很少）

Ifuel = 602.07（8.9的中性轴转动惯量占比很少）

注意driver和fuel的转动惯量实际上几乎由平行轴决定，组件的中性轴占比很少

转动惯量是给转动服务的，转动需要一个量去描述它，我们用开普勒第二定律和牛顿定律一起，类比动量的概念，推广出角动量 H = r*m*v，用角速度ω替换速度，得到H = m*r x ω x r，用转动惯性I替换，得到H = ω*I，对时间进行微分后，可以得到角动量随时间的变化率M = I*α （α是物体相对于给定轴的角加速度），角动量的随时间变化率，实际上也就是物体力矩的总和（扭矩）

在这条公式里，角动量变化率M和角加速度α都是向量，但转动惯量I会随转动轴的改变而改变，是张量，张量是指一种具有大小和方向的数学表达式，大小依据方向且不唯一。我们下面会得知I是二阶张量，3x3的矩阵

但在这之前，我们要先处理一个坐标系的问题。我们在上面计算的所有数值都是在世界坐标系下完成的，这导致在物体位置改变时，转动惯量也会发生改变。

这个可以参考这篇报告https://wenku.baidu.com/view/e8f0b36124c52cc58bd63186bceb19e8b8f6ecc4.html，无论物体在世界坐标系是发生了平移还是旋转，转动惯量都会发生改变

为了减少计算，我们需要得到在局部坐标系下的公式。我这里直接给出了，有兴趣的可以自己百度

H = I(dω / dt) + (ω x (Iω))

解决了坐标系的问题，我们可以开始求解转动惯量了，我再重申一次，转动惯量是随着转动轴变化而变化大小的张量。我们上面解决的是坐标系不同导致转动惯量的计算矩阵不同的问题。下面要解决的，是如何得到一个3x3矩阵，代表在给定的转动轴下，该矩阵乘以角速度可以得到角动量

熟悉线性代数的都知道，三维空间下的任意轴的转动都可以用一个3x3的矩阵来表示，所以很显然的想法，计算矩阵也应该是一个3x3的矩阵

如果你连3x3矩阵这点都存在疑问，也没关系，我们只需要把角动量H = m * r x ω x r展开写出来，自然可以得到以下矩阵I，使得满足以下等式

其中惯性积（Ixx,Iyy,Izz叫做转动惯量，其他的叫做惯性积）的值如下：

Ixy = Iyx = ∫(xy)dm

Ixz = Izx = ∫(xz)dm

Iyz = Izy = ∫(yz)dm

这就是转动惯量的计算矩阵了。

在有些教材里，会默认只计算只有主轴的简单物体（物体可以沿着某一坐标轴对称，则称这条轴为主轴），这些物体的惯性积都是0。但是我们在做游戏的时候，很有可能不会使用物体的主轴，所以还是有必要了解下惯性积是什么东西。

惯性积有类似于平行轴一样的公式

Ixy = I0xy + mdxdy

I0xy表示局部惯性积，即相对于通过物体重心的轴的惯性积，如果你的元素是简单的几何体，那么局部惯性积会是0（这也就是为什么我们总是用简单几何体去近似组件）

最后我附上计算转动惯量的计算矩阵的代码，方便大家理解上述过程。注意这段代码是相对于通过组合刚体重心的轴的，所以应当预先处理好组件的本地坐标，让其进入到相对于组合刚体重心的坐标系

typedef struct _PointMass
{
	float mass;
	Vector designPosition;  //全局坐标
	Vector correctedPosition; //本地坐标
	Vector localIntertia; //局部转动惯量
} PointMass;

float Ixx, Iyy, Izz, Ixy, Ixz, Iyz;
Matrix3x3 InertiaTensor;

Ixx = 0; Iyy = 0; Izz = 0;
Ixy = 0; Ixz = 0; Iyz = 0;

for (i = 0; i < _NUMELEMENTS; i++) {
	Ixx += Element[i].localIntertia.x +
		Element[i].mass * (Element[i].correctedPosition.y * Element[i].correctedPosition.y + 
		Element[i].correctedPosition.z * Element[i].correctedPosition.z); 
		
	Iyy += Element[i].localIntertia.y +
		Element[i].mass * (Element[i].correctedPosition.z * Element[i].correctedPosition.z + 
		Element[i].correctedPosition.x * Element[i].correctedPosition.x); 
		
	Izz += Element[i].localIntertia.z +
		Element[i].mass * (Element[i].correctedPosition.x * Element[i].correctedPosition.x + 
		Element[i].correctedPosition.y * Element[i].correctedPosition.y); 
		
	Ixy += 
		Element[i].mass * (Element[i].correctedPosition.x * Element[i].correctedPosition.y); 
		
	Ixz += 
		Element[i].mass * (Element[i].correctedPosition.x * Element[i].correctedPosition.z);

	Iyz += 
		Element[i].mass * (Element[i].correctedPosition.y * Element[i].correctedPosition.z); 
}

InertiaTensor.e11 =  Ixx;
InertiaTensor.e12 = -Ixy;
InertiaTensor.e13 = -Ixz;
InertiaTensor.e21 = -Ixy;
InertiaTensor.e22 = Iyy;
InertiaTensor.e23 = -Iyz;
InertiaTensor.e31 = -Ixz;
InertiaTensor.e32 = -Iyz;
InertiaTensor.e33 = Izz;

非恒定加速度运动

一般来说，如果我们需要模拟阻力，则需要解决非恒定加速度的计算问题。阻力一般与速度平方成正比，质量不变的的话，阻力加速度也与速度平方成正比。我下面给出计算式的推导过程，不感兴趣的可以往下翻结论

如果S1=0，即物体初始位移为0，则S = ln(1+vkt)/t

粘滞阻力

理想的粘滞阻力是速度和一个测出来的阻力系数的函数，阻力系数受表面情况，流体属性，流体流动状况决定，一般可以表示为： F = Cf * v，Cf是阻力系数，v是物体速度。这个公式适合在液体中低速移动的物体，这说明物体周围的流体流动是层状的，流体的流线是平行的，未被打乱。对于快速的移动的物体，粘滞阻力是：F = Cf * v^2，这说明物体周围的流体被扰动了。空气也是流体

3D Game Programming Design(六)：物理系统与碰撞

猜你喜欢