从零实现一个2D刚体物理引擎：数值积分、SAT碰撞检测与Sequential Impulse约束求解
在上一篇文章中，我们从零实现了一个Markdown解析器，完成了块级状态机和行内分隔符栈的完整链路。那篇文章聚焦的是文本结构的解析。这一次我们把视角转向数值计算和几何算法，目标是一个更“物理”的系统：2D刚体物理引擎。
物理引擎是游戏引擎、机器人仿真、交互式动画的底层支撑。它看起来复杂，但核心循环可以拆解为三个相对独立的问题：积分（给定力和速度，求解下一时刻的位置）、碰撞检测（判断两个物体是否相交，以及相交的具体信息）、约束求解（消除穿透、施加正确的冲量来模拟碰撞响应）。matter.js、p2.js、planck.js这些成熟的JS物理引擎，本质上都是在解决这三个问题，差异在于算法的精度、稳定性和工程复杂度。
本文用纯前端JavaScript从零实现一个支持凸多边形和圆的2D刚体物理引擎，不依赖任何库。完整链路是：Semi-implicit Euler积分 → SAT分离轴碰撞检测 → 接触流形生成 → Sequential Impulse约束求解 → Canvas渲染。代码可直接在浏览器中运行。
一、物理引擎的核心循环
物理引擎的主循环可以概括为“检测→求解→积分”三步。每一步的输入输出如下：
检测阶段：给定当前所有物体的位置和旋转，找出所有相互接触的物体对，并为每一对生成接触流形（contact manifold）。接触流形包含接触点的位置、穿透深度和接触法线。
求解阶段：给定接触流形和物体的质量、惯量、速度，计算需要施加的冲量，使得物体在接触点处不再相互穿透，同时满足摩擦和恢复系数约束。
积分阶段：根据求解得到的冲量更新物体速度，然后根据速度更新位置和旋转。同时施加外力（如重力）并做阻尼处理。
planck.js的核心概念文档中把World描述为“bodies、fixtures和constraints的集合”，Solver负责“推进时间并求解接触和关节约束”。这个描述准确概括了物理引擎的职责边界：World是状态容器，Solver是算法核心，积分是状态推进。
二、刚体的状态表示与数值积分
一个2D刚体在平面上的状态可以用位置、旋转角、线速度和角速度来完整描述。质量属性包括质量、质心位置和转动惯量。
javascript
class RigidBody {
constructor(options = {}) {
// 位置与旋转
this.position = { x: 0, y: 0 };
this.angle = 0;
this.velocity = { x: 0, y: 0 };
this.angularVelocity = 0;
this.force = { x: 0, y: 0 };
this.torque = 0;
// 质量属性
this.mass = options.mass ?? 0;
this.invMass = this.mass > 0 ? 1 / this.mass : 0;
this.inertia = options.inertia ?? 0;
this.invInertia = this.inertia > 0 ? 1 / this.inertia : 0;
// 形状
this.shape = options.shape ?? null;
this.restitution = options.restitution ?? 0.2; // 恢复系数
this.friction = options.friction ?? 0.3; // 摩擦系数
this.isStatic = options.isStatic ?? false;
}
applyForce(fx, fy) {
this.force.x += fx;
this.force.y += fy;
}
applyImpulse(ix, iy, contactVector = { x: 0, y: 0 }) {
// 线速度变化
this.velocity.x += ix * this.invMass;
this.velocity.y += iy * this.invMass;
// 角速度变化：冲量 × 接触点相对质心的叉积
this.angularVelocity += this.invInertia *
(contactVector.x * iy - contactVector.y * ix);
}
}
积分使用Semi-implicit Euler方法：先用当前加速度更新速度，再用更新后的速度更新位置。与显式Euler相比，Semi-implicit Euler在弹簧和碰撞场景下更稳定，因为它引入了数值阻尼，能抑制能量的指数增长。
javascript
function integrate(body, dt, gravity) {
if (body.isStatic || body.invMass === 0) return;
// 施加外力（重力）
body.force.x += gravity.x * body.mass;
body.force.y += gravity.y * body.mass;
// 1. 用当前加速度更新速度
body.velocity.x += body.force.x * body.invMass * dt;
body.velocity.y += body.force.y * body.invMass * dt;
body.angularVelocity += body.torque * body.invInertia * dt;
// 2. 阻尼（模拟空气阻力，防止速度无限增长）
body.velocity.x *= 0.998;
body.velocity.y *= 0.998;
body.angularVelocity *= 0.998;
// 3. 用更新后的速度更新位置
body.position.x += body.velocity.x * dt;
body.position.y += body.velocity.y * dt;
body.angle += body.angularVelocity * dt;
// 4. 清除累积的力和力矩
body.force.x = 0;
body.force.y = 0;
body.torque = 0;
}
转动惯量取决于形状。对于质量为m、边长为s的正方形，绕质心的转动惯量为m * s^2 / 6。对于半径为r的圆，为m * r^2 / 2。这些值在创建刚体时根据形状自动计算。
三、碰撞检测：SAT分离轴定理
碰撞检测分为宽相和窄相两个阶段。宽相用AABB（轴对齐包围盒）快速排除明显不相交的物体对，窄相用精确算法判断是否真正相交并生成接触信息。本文实现的是窄相中的核心算法——分离轴定理（Separating Axis Theorem，SAT）。
SAT的数学基础很简单：两个凸多边形不相交，当且仅当存在一条轴，使得两个多边形在该轴上的投影区间不重叠。对于凸多边形，候选轴的数量是有限的——只需要检查每个多边形的每条边的法线方向。
javascript
function getVertices(body) {
// 返回世界坐标系下的顶点列表
const cos = Math.cos(body.angle);
const sin = Math.sin(body.angle);
return body.shape.vertices.map(v => ({
x: body.position.x + v.x * cos - v.y * sin,
y: body.position.y + v.x * sin + v.y * cos
}));
}
function projectVertices(vertices, axis) {
let min = Infinity, max = -Infinity;
for (const v of vertices) {
const proj = v.x * axis.x + v.y * axis.y;
if (proj < min) min = proj;
if (proj > max) max = proj;
}
return { min, max };
}
function satPolygons(polyA, polyB) {
const vertsA = getVertices(polyA);
const vertsB = getVertices(polyB);
const axes = [];
// 收集所有边的法线作为候选轴
for (let i = 0; i < vertsA.length; i++) {
const j = (i + 1) % vertsA.length;
const edge = { x: vertsA[j].x - vertsA[i].x, y: vertsA[j].y - vertsA[i].y };
axes.push({ x: -edge.y, y: edge.x }); // 法线
}
for (let i = 0; i < vertsB.length; i++) {
const j = (i + 1) % vertsB.length;
const edge = { x: vertsB[j].x - vertsB[i].x, y: vertsB[j].y - vertsB[i].y };
axes.push({ x: -edge.y, y: edge.x });
}
let minOverlap = Infinity;
let bestAxis = null;
for (const axis of axes) {
// 归一化轴向量
const len = Math.sqrt(axis.x * axis.x + axis.y * axis.y);
const normAxis = { x: axis.x / len, y: axis.y / len };
const projA = projectVertices(vertsA, normAxis);
const projB = projectVertices(vertsB, normAxis);
// 检查投影是否重叠
if (projA.max < projB.min || projB.max < projA.min) {
return null; // 存在分离轴，不相交
}
// 计算重叠深度
const overlap = Math.min(projA.max - projB.min, projB.max - projA.min);
if (overlap < minOverlap) {
minOverlap = overlap;
bestAxis = normAxis;
}
}
// 确定法线方向：从A指向B
const centerA = polyA.position;
const centerB = polyB.position;
const d = { x: centerB.x - centerA.x, y: centerB.y - centerA.y };
if (d.x * bestAxis.x + d.y * bestAxis.y < 0) {
bestAxis.x = -bestAxis.x;
bestAxis.y = -bestAxis.y;
}
return { normal: bestAxis, depth: minOverlap };
}
SAT返回的depth是最小平移距离——如果要把两个多边形分开，沿着normal方向平移depth距离即可。这个信息在约束求解中直接用于计算恢复冲量。
圆与凸多边形的碰撞检测需要特殊处理。圆不是多边形，没有边。做法是：找到多边形上距离圆心最近的顶点或边，计算圆心到该点的距离。如果距离小于半径，则发生碰撞。
javascript
function satCirclePolygon(circleBody, polyBody) {
const center = circleBody.position;
const radius = circleBody.shape.radius;
const verts = getVertices(polyBody);
// 找到多边形上距离圆心最近的边
let minDist = Infinity;
let closestPoint = null;
for (let i = 0; i < verts.length; i++) {
const j = (i + 1) % verts.length;
const point = closestPointOnSegment(center, verts[i], verts[j]);
const dist = distance(center, point);
if (dist < minDist) {
minDist = dist;
closestPoint = point;
}
}
if (minDist >= radius) return null;
const normal = normalize({
x: closestPoint.x - center.x,
y: closestPoint.y - center.y
});
return {
normal,
depth: radius - minDist,
contactPoint: closestPoint
};
}
四、接触流形与Sequential Impulse求解器
SAT给出的是几何信息——法线和深度。约束求解需要的是接触流形：在接触点处，两个物体的相对速度、有效质量，以及需要施加的冲量。
Sequential Impulse是Erin Catto（Box2D作者）提出的迭代求解方法，本质上是投影高斯-赛德尔（Projected Gauss-Seidel） 的一种矩阵自由形式实现。Bullet物理引擎的btSequentialImpulseConstraintSolver正是采用了这一方法，官方文档描述它为“Projected Gauss Seidel (iterative LCP) method”的快速实现。
算法的核心思想是：对每个接触点，计算一个冲量，使得接触点处的相对法向速度变为零（或按恢复系数反弹）。然后迭代这个过程多次，每次迭代都会修正之前冲量引入的误差，最终收敛到满足所有接触约束的速度场。
javascript
function solveContact(contact, dt) {
const { bodyA, bodyB, normal, depth, contactPoint } = contact;
// 接触点相对质心的向量
const ra = { x: contactPoint.x - bodyA.position.x, y: contactPoint.y - bodyA.position.y };
const rb = { x: contactPoint.x - bodyB.position.x, y: contactPoint.y - bodyB.position.y };
// 接触点处的相对速度
const relativeVelocity = {
x: (bodyB.velocity.x - bodyB.angularVelocity * rb.y) -
(bodyA.velocity.x - bodyA.angularVelocity * ra.y),
y: (bodyB.velocity.y + bodyB.angularVelocity * rb.x) -
(bodyA.velocity.y + bodyA.angularVelocity * ra.x)
};
// 法向相对速度
const velAlongNormal = relativeVelocity.x * normal.x + relativeVelocity.y * normal.y;
// 如果物体正在分离，不施加冲量
if (velAlongNormal > 0) return;
// 恢复系数
const e = Math.min(bodyA.restitution, bodyB.restitution);
// 计算有效质量（标量）
const raCrossN = ra.x * normal.y - ra.y * normal.x;
const rbCrossN = rb.x * normal.y - rb.y * normal.x;
const invMassSum = bodyA.invMass + bodyB.invMass +
raCrossN * raCrossN * bodyA.invInertia +
rbCrossN * rbCrossN * bodyB.invInertia;
// 施加的冲量大小
const j = -(1 + e) * velAlongNormal / invMassSum;
// 应用冲量
applyImpulse(bodyA, -j * normal.x, -j * normal.y, ra);
applyImpulse(bodyB, j * normal.x, j * normal.y, rb);
// 位置修正（Baumgarte稳定化）：消除穿透
const percent = 0.4; // 修正比例
const slop = 0.01; // 允许的微小穿透
const correction = Math.max(depth - slop, 0) / invMassSum * percent;
bodyA.position.x -= correction * normal.x * bodyA.invMass;
bodyA.position.y -= correction * normal.y * bodyA.invMass;
bodyB.position.x += correction * normal.x * bodyB.invMass;
bodyB.position.y += correction * normal.y * bodyB.invMass;
}
function solveContacts(contacts, dt, iterations = 10) {
for (let i = 0; i < iterations; i++) {
for (const contact of contacts) {
solveContact(contact, dt);
}
}
}
这段代码中有几个关键设计。velAlongNormal > 0的判断避免了在物体正在分离时施加不必要的冲量，这是维持稳定性的重要条件。invMassSum的计算包含了角惯量的贡献——当冲量施加在偏离质心的位置时，一部分能量会转化为角速度，有效质量因此增大，所需的冲量变小。
Baumgarte稳定化是消除穿透的标准技术。纯速度求解无法完全消除穿透，因为速度冲量只能改变未来的运动，不能修正已经发生的重叠。correction项直接将物体沿法线方向推开一个与穿透深度成正比的量，percent控制修正强度（通常0.2到0.8），slop允许微小的穿透以避免抖动。
五、摩擦力的处理
摩擦力是物理引擎中最难处理的部分之一。标准的做法是在求解法向冲量之后，沿接触面的切线方向施加一个冲量，其大小受库仑摩擦定律约束：|j_t| <= μ * |j_n|。
javascript
function solveFriction(contact, dt) {
const { bodyA, bodyB, normal, contactPoint } = contact;
const ra = { x: contactPoint.x - bodyA.position.x, y: contactPoint.y - bodyA.position.y };
const rb = { x: contactPoint.x - bodyB.position.x, y: contactPoint.y - bodyB.position.y };
// 计算切线方向（法线逆时针旋转90度）
const tangent = { x: -normal.y, y: normal.x };
const relativeVelocity = {
x: (bodyB.velocity.x - bodyB.angularVelocity * rb.y) -
(bodyA.velocity.x - bodyA.angularVelocity * ra.y),
y: (bodyB.velocity.y + bodyB.angularVelocity * rb.x) -
(bodyA.velocity.y + bodyA.angularVelocity * ra.x)
};
const velAlongTangent = relativeVelocity.x * tangent.x + relativeVelocity.y * tangent.y;
const raCrossT = ra.x * tangent.y - ra.y * tangent.x;
const rbCrossT = rb.x * tangent.y - rb.y * tangent.x;
const invMassSum = bodyA.invMass + bodyB.invMass +
raCrossT * raCrossT * bodyA.invInertia +
rbCrossT * rbCrossT * bodyB.invInertia;
let jt = -velAlongTangent / invMassSum;
// 库仑摩擦约束：切向冲量不超过法向冲量 × 摩擦系数
const mu = Math.sqrt(bodyA.friction * bodyB.friction);
// 这里需要知道上一次法向冲量的大小，简化处理用深度近似
const maxFriction = mu * Math.abs(contact.normalImpulse || 1);
jt = Math.max(-maxFriction, Math.min(maxFriction, jt));
applyImpulse(bodyA, -jt * tangent.x, -jt * tangent.y, ra);
applyImpulse(bodyB, jt * tangent.x, jt * tangent.y, rb);
}
工程实现中，摩擦冲量的上限需要参考法向冲量。实际引擎（如Box2D、planck.js）会在求解器中存储每个接触点的法向冲量累积值，摩擦求解时引用这个值。上面的简化实现用contact.normalImpulse近似，在完整实现中需要在solveContact中把j存入contact对象。
六、碰撞检测的宽相优化
上面的SAT是窄相算法，对每一对物体都执行一次。如果场景中有N个物体，朴素做法是O(N²)对检测。对于几十个物体这还可以接受，但对于几百个物体的场景，宽相优化是必需的。
最简单的宽相优化是Sweep and Prune（排序扫描）算法。它的核心思想是：如果两个物体在x轴上的投影区间不重叠，它们一定不相交。将所有物体按AABB的最小x坐标排序，然后只需要检查排序后相邻的物体对。
javascript
function broadPhase(bodies) {
// 1. 计算每个物体的AABB
const entries = bodies.map((body, index) => {
const aabb = computeAABB(body);
return { body, index, minX: aabb.minX, maxX: aabb.maxX };
});
// 2. 按minX排序
entries.sort((a, b) => a.minX - b.minX);
// 3. 扫描：对每个物体，只检查minX小于其maxX的后续物体
const pairs = [];
for (let i = 0; i < entries.length; i++) {
const a = entries[i];
for (let j = i + 1; j < entries.length; j++) {
const b = entries[j];
if (b.minX > a.maxX) break; // 后续的minX更大，不可能重叠
// 检查y轴重叠
const aabbA = computeAABB(a.body);
const aabbB = computeAABB(b.body);
if (aabbA.minY <= aabbB.maxY && aabbB.minY <= aabbA.maxY) {
pairs.push([a.body, b.body]);
}
}
}
return pairs;
}
Sweep and Prune将宽相复杂度从O(N²)降低到O(N log N + K)，其中K是实际重叠的物体对数。在物体稀疏分布的典型场景中，K远小于N²。
p2.js作为成熟的JS物理引擎，其功能包括“collision detection, contacts, friction, restitution, motors, springs, advanced constraints and various shape types”。我们的实现覆盖了其中碰撞检测、接触、摩擦和恢复系数的核心部分，约束类型只实现了接触约束，还没有实现距离约束、旋转约束等关节类型。
七、完整物理世界与渲染
把所有部分串联起来，物理世界的更新循环如下：
javascript
class PhysicsWorld {
constructor() {
this.bodies = [];
this.gravity = { x: 0, y: 9.81 * 100 }; // 像素尺度下的重力
this.contacts = [];
}
addBody(body) { this.bodies.push(body); }
step(dt) {
// 1. 宽相检测
const pairs = broadPhase(this.bodies.filter(b => !b.isStatic));
// 2. 窄相检测 + 接触流形生成
this.contacts = [];
for (const [a, b] of pairs) {
const contact = detectCollision(a, b);
if (contact) this.contacts.push(contact);
}
// 3. 约束求解
solveContacts(this.contacts, dt, 10);
// 4. 积分
for (const body of this.bodies) {
integrate(body, dt, this.gravity);
}
}
}
渲染用Canvas 2D API。多边形用路径描边和填充，圆用arc。为了可视化旋转，可以在多边形上画一条从质心到某个顶点的线段。
javascript
function render(ctx, world) {
ctx.clearRect(0, 0, ctx.canvas.width, ctx.canvas.height);
for (const body of world.bodies) {
ctx.save();
ctx.translate(body.position.x, body.position.y);
ctx.rotate(body.angle);
if (body.shape.type === 'circle') {
ctx.beginPath();
ctx.arc(0, 0, body.shape.radius, 0, Math.PI * 2);
} else {
ctx.beginPath();
const verts = body.shape.vertices;
ctx.moveTo(verts[0].x, verts[0].y);
for (let i = 1; i < verts.length; i++) {
ctx.lineTo(verts[i].x, verts[i].y);
}
ctx.closePath();
}
ctx.fillStyle = body.color || '#4a90d9';
ctx.fill();
ctx.strokeStyle = '#222';
ctx.stroke();
// 绘制质心标记
ctx.beginPath();
ctx.arc(0, 0, 3, 0, Math.PI * 2);
ctx.fillStyle = '#e74c3c';
ctx.fill();
ctx.restore();
}
}
创建一个简单的场景来测试：
javascript
const world = new PhysicsWorld();
// 地面（静态）
world.addBody(new RigidBody({
shape: { type: 'polygon', vertices: [
{ x: -400, y: 300 }, { x: 400, y: 300 },
{ x: 400, y: 350 }, { x: -400, y: 350 }
]},
isStatic: true
}));
// 几个掉落的方块和圆
for (let i = 0; i < 5; i++) {
const size = 40;
world.addBody(new RigidBody({
shape: { type: 'polygon', vertices: [
{ x: -size/2, y: -size/2 }, { x: size/2, y: -size/2 },
{ x: size/2, y: size/2 }, { x: -size/2, y: size/2 }
]},
mass: 1,
inertia: 1 * size * size / 6,
position: { x: -100 + i * 60, y: -200 },
color: `hsl(${i * 60}, 70%, 60%)`
}));
}
// 主循环
let lastTime = performance.now();
function loop(now) {
const dt = Math.min((now - lastTime) / 1000, 1/30); // 限制最大步长
lastTime = now;
world.step(dt);
render(ctx, world);
requestAnimationFrame(loop);
}
requestAnimationFrame(loop);
八、稳定性问题与工程考量
物理引擎的工程实现中，稳定性是最大的挑战。几个关键问题：
时间步长与子步。过大的dt会导致穿透和能量爆炸。标准做法是将物理更新固定为1/60秒，如果实际帧间隔更大，执行多次子步。planck.js的文档中特别提到了“sub-stepping solver”，它会“将物体移动到首次碰撞时间，然后求解碰撞”。
速度阈值与休眠。快速运动的物体会穿透薄壁（隧穿效应）。最简单的缓解是限制最大速度。更完整的方案是连续碰撞检测（CCD），在物体运动的路径上检测碰撞，而不是只检测离散时刻的位置。
堆叠稳定性。当多个物体叠在一起时，每个接触点的求解都会影响相邻的接触点，导致“抖动”。Sequential Impulse的迭代次数不够时，堆叠会缓慢下沉或抖动。增加迭代次数到15到20可以显著改善，但代价是性能。工程上通常用10次迭代，配合Baumgarte稳定化和位置修正来缓解。
角惯量的计算。对于任意凸多边形，转动惯量需要根据顶点坐标积分计算，而不是简单的公式。标准做法是将多边形分解为三角形，然后累加每个三角形对质心的惯量贡献。
九、总结
从刚体状态表示到Semi-implicit Euler积分，从SAT分离轴碰撞检测到接触流形生成，从Sequential Impulse约束求解到Baumgarte位置修正——这个2D物理引擎的核心代码不到400行，但覆盖了实时物理仿真的全部关键算法。
SAT的数学基础是凸集的分离定理，Sequential Impulse是投影高斯-赛德尔方法的矩阵自由实现，Semi-implicit Euler是辛积分器的一种。这些算法在游戏引擎、机器人仿真、交互式动画中反复出现，理解了它们的数学原理和工程实现，再去阅读Box2D、Bullet或PhysX的源码，会发现核心逻辑是相通的，差异只在优化策略——SIMD并行化、约束图的岛屿分割、连续碰撞检测的TOI算法。物理引擎的工程复杂度远远超出算法本身，但算法是所有工程优化的起点。
- 点赞
- 收藏
- 关注作者
评论（0）