基于动量计算的刚体碰撞模拟【Games103 |
您所在的位置：网站首页 › 刚体动量如何计算 › 基于动量计算的刚体碰撞模拟【Games103 |
基于动量计算的刚体碰撞模拟【Games103
|
Games103的lab1的学习笔记和代码实现
课程链接： GAMES103-基于物理的计算机动画入门_哔哩哔哩_bilibili 课程主页：GAMES103：基于物理的计算机动画入门 – 计算机图形学与混合现实在线平台 1. 理论基础 1.1 点乘与叉乘 1.2 动量与角动量转动惯量是刚体绕轴转动时惯性的量度，可以理解为线性动力学中的质量，是一个标量，不需要转换空间。
为了避免万向节死锁问题（某种情况下降低旋转自由度），引入了四元数的概念。 四元数是一个超复数，其（x，y，z，w）四个值分为三个虚部，一个实部。 与欧拉角的转换可以参考这个文章3D数学基础:四元数与欧拉角之间的转换 - rainbow70626 - 博客园
即xyz三个值是在xyz轴的分量，确定了一个旋转轴的方向，然后w是绕旋转轴旋转的角度。 1.5 坐标空间转换 2. 算法实现 2.1 算法思路本质上刚体发生变化的就是位置（向量）和旋转（四元数）。 碰撞过程导致了受力变化，从而影响物体速度，进而更新物体的位置坐标和旋转。本算法是当前帧如果发生碰撞，计算出物体局部受到冲量大小（没有发生但是假设发生），从而求出整个物体的速度变化，然后对位置更新（更新当前帧）。优点是可以及时对受力进行反馈。
所以如果发生碰撞，需要对于局部受碰撞点，受到弹力和摩擦力作用，计算出之前的和新的速度大小，得到局部冲量大小。 2.2 刚体碰撞 2.2.1 碰撞检测检测方法是非常直观的，计算点到表面的距离，大于则在外面，小于在里面。难点是对于复杂网格表面，应该如何计算。
遍历每个点（世界坐标），计算到碰撞体表面的距离（PXi与法线点乘），负数表示碰撞。
对于局部受碰撞点： 计算法向和切向速度，法向受弹力作用：速度大小近似表示为法线的反向乘弹力系数 切向受到摩擦力作用：按照a系数减少，但是不能小于0（速度不会反方向回去） 课程表示a系数需要满足库伦定律，具体a的计算如下图所示。
对于受碰撞的物体： 如果我们知道冲量大小，那么就可以求到新的速度和角速度的值。
所以我们需要对局部碰撞点求得冲量大小： 得到计算式如下：
星号表示向量的叉乘可以用一个新的矩阵表示，简化了很多（tql）
知道速度和角速度后我们就可以很快求出新的位置和旋转。
1. 碰撞点不止一个，可能有几百上千个，所以我们求的碰撞点是它们的平均坐标 2. 碰撞时候会因为一直有向下的重力加速度，所以弹力也一直存在，导致抖动问题，因此在趋于稳定的时候，将弹力减小，使得速度尽可能控制为0 3. unity里面没有四元数的加法，unity中的四元数表示如下： 具体代码 using UnityEngine; using System.Collections; public class Rigid_Bunny : MonoBehaviour { bool launched = false; float dt = 0.015f; Vector3 v = new Vector3(0, 0, 0); // velocity Vector3 w = new Vector3(0, 0, 0); // angular velocity Vector3 gravity = new Vector3(0,-9.8f,0); float mass; // mass Matrix4x4 I_ref; // reference inertia float linear_decay = 0.999f; // for velocity decay float angular_decay = 0.98f; float restitution = 0.3f; // for collision float friction = 0.2f; // 摩擦系数 bool lowCollide = false; Mesh mesh; // Use this for initialization void Start () { transform.rotation = new Quaternion(0, Mathf.Sin(Mathf.Deg2Rad * 30) * 1, 0,Mathf.Sin(Mathf.Deg2Rad * 30)); mesh = GetComponent().mesh; Vector3[] vertices = mesh.vertices; float m=1; mass=0; for (int i=0; i 0) { // campute the wanted vi_new avgPos /= validNum; Vector3 r_i = avgPos; Vector3 v_i = v + Vector3.Cross(w, R.MultiplyVector(r_i)); if(Vector3.Dot(v_i,N) > 0)//the velocity is positive return; Vector3 x_i = x + R.MultiplyVector(r_i); float distance = Vector3.Dot(x_i, N); print("distance" + distance + "posY" + x_i.y); if (-0.3f < distance) { lowCollide = true; } Vector3 vi_normal = Vector3.Dot(v_i,N) * N; Vector3 vi_tangent = v_i - vi_normal; Vector3 vi_normal_new = -restitution * vi_normal; float a = Mathf.Max(1- friction * (1+restitution)*(vi_normal.magnitude/vi_tangent.magnitude),0); Vector3 vi_tangent_new = a * vi_tangent; Vector3 vi_new = vi_normal_new + vi_tangent_new; // compute the impulse //Matrix4x4 I_world = R * I_ref * Matrix4x4.Transpose(R);// inertia is a quantity Matrix4x4 I_inv = I_ref.inverse; Matrix4x4 Rc = Get_Cross_Matrix(R * r_i); Matrix4x4 K = Matrix_Subtract(MatrixMultiplyFloat(Matrix4x4.identity,1/mass) , Rc * I_inv * Rc); Vector3 j = K.inverse * (vi_new - v_i); // update v and w Vector3 dv = j / mass; Debug.DrawRay(x + R.MultiplyVector(r_i), dv, Color.red);// draw the velocity direction //Debug.DrawRay(x,v,Color.green);// gobal velocity v += dv; w = w + (Vector3)(I_inv * Vector3.Cross(R * r_i , j)); } } Quaternion QuaternionAdd(Quaternion q1,Quaternion q2) { q1.x += q2.x; q1.y += q2.y; q1.z += q2.z; q1.w += q2.w; return q1; } // Update is called once per frame void Update () { //Game Control if(Input.GetKey(KeyCode.R)) { transform.position = new Vector3 (0, 2f, 0); restitution = 0.5f; launched=false; } if(Input.GetKey(KeyCode.S))// shoot { v = new Vector3 (5, 2, 0); launched=true; } // Part I: Update velocities v += gravity * dt; v *= linear_decay; w *= angular_decay; // Part II: Collision Impulse lowCollide = false; Collision_Impulse(new Vector3(0, 0.01f, 0), new Vector3(0, 1, 0)); Collision_Impulse(new Vector3(2, 0, 0), new Vector3(-1, 0, 0)); if (lowCollide && v.magnitude < 0.13f) { restitution *= 0.5f; if (restitution < 0.01f) { restitution = 0; } if(v.magnitude < 0.001f) { v = Vector3.zero; } } else { restitution = 0.5f; } v *= linear_decay; w *= angular_decay; // Part III: Update position & orientation //Update linear status Vector3 x = transform.position; x += v * dt; //Update angular status Quaternion q = transform.rotation; Quaternion qw = new Quaternion(w.x * dt/2, w.y * dt / 2, w.z * dt / 2, 0); q = QuaternionAdd(q , qw * q); // Part IV: Assign to the object transform.position = x; transform.rotation = q; } } |
【本文地址】
今日新闻 |
推荐新闻 |