物理引擎中可能是浮点数误差造成的一个问题?

发布于 2022-09-01 05:15:54 字数 3629 浏览 24 评论 0

自己试着做个2d物理引擎玩时遇到的一个问题... 两个物体碰撞时计算摩擦冲量似乎很容易受浮点数精度的影响

两个物体碰撞时给予对方的碰撞冲量我按下面这个公式算了出来

其中速度差向量 vr = vA' - vB' , n是两物体分离轴的法线单位向量
另外把角速度也换算成线速度加到速度向量上来计算, 于是vA' = vA + rA×ωA
到这里都没什么问题

现在要求摩擦冲量, 令单位切线向量

$$tangent = normalize(vr - (vr*normal) * normal)$$

于是变成下图这样

按图1的上面的公式求得碰撞冲量大小j之后, 摩擦冲量大小$$jf = μ * j$$, 摩擦力在两个物体相对静止时就会变成0所以还要求一个最大摩擦冲量大小
vr在tangent上的投影的大小 $$vt = vr * tangent$$, 对两个物体施加摩擦冲量之后最多使得这个vt变成0但不能让vt的符号改变

我按照图1下面的4个公式以及△v = vt求得的最大jf是:
$$jfMax = vr * tangent / (1/mA + 1/mB + (rA×tangent)^2/IA + (rB×tangent)^2/IB)$$
于是摩擦冲量 $$tImpulse = tangent * min(jfMax, μ * j)$$

分别对物体A施加 -tImpulse, 对物体B施加 tImpulse, 如此实现摩擦

但是我发现在去除μ*j只适用jfMax的情况下, 最后vt并不等于0而是在0附近上下波动, 多数是0.01~0.5

bodyA.applyImpulse(tangentImpulse.negate(), cp);
bodyB.applyImpulse(tangentImpulse, cp);
vA = bodyA.velocity.add(rA.scalarCross(bodyA.angularVelocity));
vB = bodyB.velocity.add(rB.scalarCross(bodyB.angularVelocity));
console.log(vA.sub(vB).dot(tangent));

适用摩擦冲量前后, 有时动能反而会增大

k1 = bodyA.getKinetic() + bodyB.getKinetic();
bodyA.applyImpulse(tangentImpulse.negate(), cp);
bodyB.applyImpulse(tangentImpulse, cp);
k2 = bodyA.getKinetic() + bodyB.getKinetic();
if ((k2 - k1) > 0.1) console.log(((k2 - k1) / k1 * 100 | 0) + '%');

单个小物体会因为这个问题在与其它物体碰撞时明显的"抽搐"

求摩擦冲量的式子我翻了下别人的物理引擎基本都一样, 所以我觉得可能是浮点数的精度误差累积造成的问题
别人的代码里有看到一些误差修正的式子, 但搞不懂是怎么得出来的也不知道实际修正的是个啥...

另外我的一大坨代码如下

applyImpulse: function () {
            var bodyA = this.bodyA,
                bodyB = this.bodyB,
                cp = this.contactPoint,
                restitution = bodyA.restitution > bodyB.restitution ? bodyA.restitution : bodyB.restitution,
                normal = this.normal,
                rA = cp.sub(bodyA.centroid),
                rB = cp.sub(bodyB.centroid),
                vA = bodyA.velocity.add(rA.scalarCross(bodyA.angularVelocity)),
                vB = bodyB.velocity.add(rB.scalarCross(bodyB.angularVelocity)),
                invIA = bodyA.inverseInertia,
                invIB = bodyB.inverseInertia,
                invMA = bodyA.inverseMass,
                invMB = bodyB.inverseMass,
                vr = vA.sub(vB),
                rsnA = rA.cross(normal),
                rsnB = rB.cross(normal),
                kn = invMA + invMB + rsnA * rsnA * invIA + rsnB * rsnB * invIB,
                j = -(1 + restitution) * vr.dot(normal) / kn,
                impulse = normal.mul(j);

            var tangent = vr.sub(normal.mul(vr.dot(normal))).normalize(),
                rstA = rA.cross(tangent),
                rstB = rB.cross(tangent),
                kt = invMA + invMB + rstA * rstA * invIA + rstB * rstB * invIB,
                jfMax = vr.dot(tangent) / kt,
                sf = Math.sqrt(bodyA.staticFriction * bodyB.staticFriction),
                df = Math.sqrt(bodyA.dynamicFriction * bodyB.dynamicFriction),
                jf = jfMax < Math.abs(j * sf) ? jfMax : Math.abs(j * df),
                tangentImpulse = tangent.mul(jf);

            bodyA.applyImpulse(impulse.add(tangentImpulse.negate()), cp);
            bodyB.applyImpulse(impulse.negate().add(tangentImpulse), cp);

            return this;
        }

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

一张白纸 2022-09-08 05:15:54

像 box2d 这类的物理引擎有明确指出在 32 位浮点范围内它只能较精準处理 10cm 到 10m 范围内的物体，超出这个范围就会出现各种诡异现象。你处理的物体的尺寸在什么范围呢？

回复收藏 0

渔村楼浪 2022-09-08 05:15:54

抱歉隔了好久才想起来... 最后发现其实只要先求出碰撞后的碰撞方向的速度, 再以此求摩擦力就对了, 浮点误差其实小得没有影响简单的乘个0.9995之类的系数就能避免摩擦力过大的, 似乎是和物体可以相互嵌入对方而实际上的刚体并不能如此有关

                var contactNum = this.contact.length,
                i;
            for (i = this.contact.length; i--;) {
                var bodyA = this.bodyA,
                    bodyB = this.bodyB,
                    cp = this.contact[i],
                    restitution = Math.sqrt(bodyA.restitution * bodyB.restitution),
                    normal = this.normal,
                    rA = cp.sub(bodyA.centroid),
                    rB = cp.sub(bodyB.centroid),
                    vA = bodyA.velocity.add(rA.scalarCross(bodyA.angularVelocity)),
                    vB = bodyB.velocity.add(rB.scalarCross(bodyB.angularVelocity)),
                    invIA = bodyA.inverseInertia,
                    invIB = bodyB.inverseInertia,
                    invMA = bodyA.inverseMass,
                    invMB = bodyB.inverseMass,
                    vr = vA.sub(vB),
                    rsnA = rA.cross(normal),
                    rsnB = rB.cross(normal),
                    kn = invMA + invMB + rsnA * rsnA * invIA + rsnB * rsnB * invIB,
                    jn = -(1 + restitution) * vr.dot(normal) / kn / contactNum,
                    impulse = normal.mul(jn);

                if(Math.abs(jn) < 0.0001) {
                    return this;
                }

                bodyA.applyImpulse(impulse, cp);
                bodyB.applyImpulse(impulse.negate(), cp);

                vA = bodyA.velocity.add(rA.scalarCross(bodyA.angularVelocity));
                vB = bodyB.velocity.add(rB.scalarCross(bodyB.angularVelocity));
                vr = vA.sub(vB);
                var tangent = vr.sub(normal.mul(vr.dot(normal))).normalize(),
                    rstA = rA.cross(tangent),
                    rstB = rB.cross(tangent),
                    kt = invMA + invMB + rstA * rstA * invIA + rstB * rstB * invIB,
                    jfMax = vr.dot(tangent) / kt / contactNum * 0.9995,
                    sf = Math.sqrt(bodyA.staticFriction * bodyB.staticFriction),
                    df = Math.sqrt(bodyA.dynamicFriction * bodyB.dynamicFriction),
                    jf = jfMax < Math.abs(jn * sf) ? jfMax : Math.abs(jn * df),
                    tangentImpulse = tangent.mul(jf);

                if(jf < 0.0001) {
                    return this;
                }

                bodyA.applyImpulse(tangentImpulse.negate(), cp);
                bodyB.applyImpulse(tangentImpulse, cp);

回复收藏 0