从零实现简易2D物理引擎：碰撞检测、冲量响应与稳定堆叠-北京尧图网络

写这篇文章之前我先把话说在前面市面上的物理引擎已经做得足够好了无论是闭源的商业引擎、开源的物理库还是游戏引擎内置的物理模块都能满足大多数场景。但你如果真的想搞懂“碰撞检测为什么这么写”“堆叠为什么总会抖”“摩擦力为什么有时候像抹了油”唯一的办法就是自己从零写一个足够简单的实现把每一行代码和背后的物理公式对上。这篇帖子就带你走完这个过程最终你会得到一个能跑、能看、能继续扩展的二维物理引擎核心。全程不依赖任何第三方图形库只要有一个能跑 Python 的环境代码就能复现。先说清楚目标我们要实现的不是 Box2D 那样的工业级物理引擎恰恰相反它足够“简易”——支持圆形和轴对齐矩形支持重力、冲量碰撞响应、库仑摩擦、基础堆叠和穿透修正。别小看这点能力市面上不少独立游戏和原型 Demo 的物理效果就是从这个量级起步的。读完并跑通代码后你会对刚体动力学、碰撞流形、迭代求解这些概念有完全不一样的体感。1. 为什么需要亲手写一个“简易物理引擎”目标、边界与引擎全景动手之前得先想明白一件事你为什么要造这个轮子。我的理由是物理引擎表面上是一个 API本质上是一套数学算法。你在使用现成引擎时遇到的大部分诡异问题——物体抖、穿透、突然弹飞——都源于对底层算法缺乏直觉。自己写一遍你会建立起一套“物理引擎是在逼近真实而不是计算精确解”的判断标准后面再去调参、优化、选型思路会清晰得多。1.1 这个引擎要解决什么问题我们针对的场景非常具体一个二维世界里面有圆球、方块、一面地面墙重力往下拉物体掉下来会互相碰撞、会堆叠、会滑动。不需要流体、软体、关节、粒子也不需要连续碰撞检测。把这些功能砍掉之后核心就只剩三件事碰撞检测判断两个形状有没有相交以及相交的深度、法线、接触点。动力学求解把重力和其它力转化为速度和角速度的变化。碰撞响应用冲量把相互穿透的物体推开同时让速度发生合理的跳变。这三件事恰好对应了物理引擎中最经典的三层结构也是很多引擎源码里模块划分的原型。从工程实现的角度看一个“简易物理引擎”的骨架大致如下面这个表模块职责本文对应章节数学库矢量、点积、叉积、归一化第2章刚体状态位置、速度、质量、转动惯量第2章动力学积分力到速度、速度到位移的推进第3章碰撞检测Broad Phase Narrow Phase第4章碰撞响应法向冲量、切向摩擦第5章稳定性处理迭代求解、穿透修正第6章这不是唯一的划分方式但作为一篇教程的骨架它足够清晰。后面的每一章都会落到这张表上的具体位置。1.2 明确“不做”的部分写代码之前还有必要把边界划清楚。我们不做旋转矩形OBB之间的通用分离轴检测只做轴对齐矩形AABB不做关节约束不做马达不做流体固定时间步长而不是自适应步长用冲量法而不是基于力的约束求解。这些简化是刻意的。轴对齐矩形的碰撞检测实现起来远比旋转矩形简单但核心思路——分离轴定理——是一样的。先跑通 AABB 版本再去扩展到 OBB路径要平滑得多。这就好比你不会一开始就给引擎加上全套优化先把一辆能跑的车造出来再谈悬挂和空气动力学。把目标订到这么具体还有一个额外的好处一旦 Demo 跑通你能亲眼看到“为什么堆叠需要迭代求解”“为什么穿透修正不能直接改位置”因为你会亲手写出那版“会塌的代码”和“修好的代码”这种对比比任何理论文章都来得直接。2. 地基工程矢量的加减乘除与刚体状态表示物理引擎里 80% 的运算发生在二维平面上力的合成、速度的分解、方向的判定。所以第一个要写好的部分就是一个可靠的二维矢量类。它不需要多高性能但要语义清晰、不犯错。我在自己的实现里一般只写六七个方法因为它们几乎构成了后面所有几何算法的最小集合。2.1 一个够用的 Vec2 实现下面这段代码是完整的二维向量类我会把它作为整个引擎的数学地基。你注意看cross方法在二维里它返回的是一个标量这个标量的大小等于 |a||b|sinθ符号表示旋转方向。这个量在冲量计算里非常关键它会反复出现。import math class Vec2: __slots__ (x, y) def __init__(self, x0.0, y0.0): self.x float(x) self.y float(y) def __add__(self, other): return Vec2(self.x other.x, self.y other.y) def __sub__(self, other): return Vec2(self.x - other.x, self.y - other.y) def __mul__(self, scalar): return Vec2(self.x * scalar, self.y * scalar) __rmul__ __mul__ def dot(self, other): return self.x * other.x self.y * other.y def cross(self, other): return self.x * other.y - self.y * other.x def length_sq(self): return self.x * self.x self.y * self.y def length(self): return math.sqrt(self.length_sq()) def normalized(self): length self.length() if length 1e-12: return Vec2(1.0, 0.0) inv 1.0 / length return Vec2(self.x * inv, self.y * inv) def perpendicular(self): # 逆时针旋转90度(-y, x) return Vec2(-self.y, self.x)为什么这么设计核心原因有三个。第一所有方法都返回新对象避免别名修改导致的隐蔽 bug第二length_sq专门用于距离比较避免无谓的开平方这在碰撞检测里是高频操作第三perpendicular方法后面摩擦和切向计算都要用它单独抽出来比每次都现写Vec2(-v.y, v.x)要清晰得多。2.2 刚体的状态字段为什么要存质量倒数有了矢量下一步是定义刚体。初学者最容易忽略的一个设计是不要直接存mass要存invMass质量的倒数。原因很实际物理引擎里大量公式形如“速度变化量冲量 × invMass”单位质量的倒数能直接把除法变成乘法效率高一点只是其次更重要的是静态物体质量无限大可以统一表示为invMass 0代码里不需要到处写if static分支。同理转动惯量的倒数invI也这样处理。转动惯量本身的计算随着形状不同而不同圆是0.5 * m * r²矩形是m * (w² h²) / 12。但在简易引擎里我建议先让圆的invI 0也就是不参与旋转这样冲量公式可以先跑通最简版本后面要开放旋转时再把对应项加回来。刚体的完整状态我用这样一组字段class Body: def __init__(self, shape, position, mass1.0, restitution0.3): self.shape shape # circle 或 aabb self.shape_data ... # 圆的半径或矩形的半宽半高 self.position position # 质心位置 self.velocity Vec2() # 线速度 self.angle 0.0 # 当前角度AABB模式下恒为0 self.angular_vel 0.0 # 角速度 self.force_acc Vec2() # 每帧重置的力累加器 self.inv_mass 0.0 if mass 0.0 else 1.0 / mass self.inv_i 0.0 # 本例中默认不旋转 self.restitution restitution # 恢复系数注意force_acc的使用习惯物理引擎里力通常是每帧清零后重新累加的这让你可以在一个帧内给物体加多个力然后在积分阶段一次性应用。如果不这么做上一帧的力会残留下来导致物体越跑越奇怪。2.3 一个容易被忽略的约定所有速度都必须相对于质心刚体的线速度是相对质心的角速度是绕质心的。在很多入门实现里这两者是分开处理的但如果你的接触点不在质心位置速度场其实应该是接触点的实际速度质心速度角速度 × (接触点 - 质心) 的垂直向量这句话现在不用完全理解只要记住一点冲量计算时如果你忽略了接触点到质心的偏移那么你算出来的相对速度就是不准确的物体会出现“明明撞上了却不转”的诡异表现。我们 Demo 里的圆形物体invI 0可以暂时跳过这个计算但概念上要留个印象。3. 运动积分半隐式欧拉为什么是默认选项物理引擎的每一帧要回答一个问题给定当前速度、加速度和时间步长物体下一秒该在哪里。这个过程的学名叫“数值积分”。你可能听说过 RK4、Verlet、各种高阶积分器的名字但在绝大多数游戏物理引擎里实际用的是最朴素的一阶方法——半隐式欧拉。3.1 显式欧拉与半隐式欧拉的差距显式欧拉是最直觉的写法位置速度 * dt 速度加速度 * dt它先更新位置再用旧速度叠加。问题是这种顺序会让系统能量漂移典型表现是弹簧越蹦越高、球越弹越远因为位置更新时用的速度已经“过时”了。半隐式欧拉把顺序反过来速度加速度 * dt 位置速度 * dt一个顺序的调换带来的稳定性提升是质的飞跃。位置更新用的是刚刚计算过的新速度相当于显式地利用了未来的信息能量漂移要小得多。你可以拿一个最简单的弹簧系统做实验分别用两种顺序跑 10 秒显式欧拉早就炸了半隐式欧拉还在规规矩矩地振荡。这也是为什么工业级引擎在默认积分器上惊人一致地选择了半隐式欧拉。积分器稳定性计算量适用场景显式欧拉差能量易扩散极低几乎不用半隐式欧拉中等足够稳定低游戏物理的默认选择Verlet中高无需存储速度低粒子系统、布娃娃RK4高高高精度科学计算简易引擎的定位决定了半隐式欧拉就是正确答案。实现同样简单得惊人def step(body, gravity, dt): # 力 - 加速度 - 速度半隐式第一步 acceleration body.force_acc * body.inv_mass body.velocity acceleration * dt body.velocity gravity * dt # 速度 - 位置半隐式第二步 body.position body.velocity * dt body.force_acc Vec2() # 清零力累加器其中gravity我一般设为Vec2(0, -9.8)朝下为负 y 方向。这样一种简单的推进方式可以看作整个物理世界的“时钟”。3.2 固定时间步长为什么不能直接每帧调一次很多初学者做引擎时会在渲染循环里直接调用step(dt)而dt是上一帧的实际耗时。这个做法有个隐蔽问题dt是波动的而物理模拟对时间步长非常敏感。同一个场景60 FPS 和 30 FPS 跑出来的结果可能完全不同因为碰撞求解对步长的依赖是线性的。正确做法是用一个累加器把渲染帧率与物理步长解耦dt_fixed 1.0 / 60.0 accumulator 0.0 def frame_update(frame_dt): global accumulator accumulator min(frame_dt, 0.25) # 防止死亡循环 while accumulator dt_fixed: physics_step(dt_fixed) accumulator - dt_fixed每次physics_step都用固定的dt_fixed模拟结果才可复现。这个技巧看起来简单但它直接影响“堆叠会不会抖、碰撞会不会穿透”的体验。建议从第一天就把这个结构搭好。3.3 积分顺序的细节一个完整物理帧的积分顺序我按下面这个顺序执行合并所有外力到force_acc用半隐式欧拉更新所有物体的速度更新所有物体的位置和角度碰撞检测生成接触对列表迭代求解所有接触对的冲量、摩擦位置修正把穿透拉回阈值内。很多人会把碰撞响应放在速度更新之前结果撞完还会穿透。顺序问题的本质在于碰撞响应处理的是“最新位置上的速度关系”所以检测应当基于最新的位置求解应当作用于最新的速度。这篇博文后面所有代码都遵循这个顺序建议你也固定下来别轻易换。4. 碰撞检测从 Broad Phase 到 SAT 分离轴原理碰撞检测是物理引擎里工作量最大、也最容易出错的部分。它分了两个阶段Broad Phase 是粗筛目标是快速排除不可能相交的物体对Narrow Phase 是精算对候选物体对做精确的几何测试输出接触流形。4.1 Broad Phase先花小代价排除大多数如果场景里的物体数量少于三四十个一种足够简单的 Broad Phase 是直接暴力遍历所有物体对每对先用包围盒做一次快速排斥测试。即使这样也能过滤掉大量不相邻的物体。包围盒用 AABB 表达就够了每个刚体维护一个min_x, max_x, min_y, max_y。两个物体要相交必须满足各自 AABB 在 x 和 y 方向都重叠。这里的优化点是AABB 可以随着物体移动每帧重算也可以惰性更新。简易引擎直接每帧重算即可因为刚体数量不多。如果想更激进一点可以用一个均匀网格或者空间哈希表把世界划分成格子每帧只检测同格子和相邻格子里的物体对。这个小优化大概几十行代码能将 n² 的复杂度降为近似 n。对于二维场景尤其有效。最怕的就是物体都挤在一个格子里退化成暴力遍历所以格子尺寸最好设置成物体平均直径的 2 倍左右。4.2 Narrow Phase三种基本形状测试Narrow Phase 是真正决定碰撞精确度的部分。我们的引擎有三个测试要写圆-圆、圆-矩形AABB、矩形-矩形AABB。下面逐个拆解它们的数学原理。圆与圆的检测最简单如果两圆心距离小于半径之和则相交。法线方向是连接圆心的方向穿透深度是半径之和减去实际距离接触点取两圆心连线的中点附近。注意不要直接开平方先用length_sq和r_sum_sq比较只在确认相交后才求平方根和法线。def circle_circle_manifold(a, b): delta b.position - a.position r_sum a.shape_data.radius b.shape_data.radius if delta.length_sq() r_sum * r_sum: return None dist delta.length() normal delta * (1.0 / dist) if dist 1e-8 else Vec2(1, 0) penetration r_sum - dist contact_point a.position normal * (a.shape_data.radius - penetration * 0.5) return normal, penetration, contact_point圆与矩形的检测思路是把圆心“夹”到矩形内部得到一个最近点然后看圆心到最近点的距离是否小于半径。这里的clamp操作是关键它把圆心映射到矩形边界上最近的位置。如果圆心在矩形内部最近点就等于圆心本身这就需要一个退化分支找出圆心离四条边最近的一面以那面的法线作为碰撞法线。def clamp_point_to_aabb(center, half_w, half_h): return Vec2( max(center.x - half_w, min(center.x, center.x)) if False else clamp(center.x, -half_w, half_w), clamp(center.y, -half_h, half_h) )为了避免把代码写拧巴我一般把 AABB 的位置字段视为“中心点”矩形的表示是center half_width half_height计算时一律先算绝对边界再与圆心比较。处理矩形内部的情况时我会显式写一个“最近面”函数把圆心相对于矩形中心的偏移分别投影到四条边的法线方向取深度最小且为正的那个作为法线。这样代码更直白也方便后续扩展旋转矩形。矩形与矩形的 AABB 测试是整个引擎里最简单也最有价值的检测。两个 AABB 相交的条件是它们在两个轴向上的重叠都大于 0。重叠较少的轴就是碰撞法线所在的轴法线方向指向重叠为正的方向。def aabb_aabb_manifold(a, b): dx abs(b.position.x - a.position.x) dy abs(b.position.y - a.position.y) ox a.shape_data.half_w b.shape_data.half_w - dx oy a.shape_data.half_h b.shape_data.half_h - dy if ox 0 or oy 0: return None if ox oy: normal Vec2(1, 0) if b.position.x a.position.x else Vec2(-1, 0) penetration ox contact_point Vec2( (a.position.x b.position.x) / 2.0, (a.position.y b.position.y) / 2.0 ) else: normal Vec2(0, 1) if b.position.y a.position.y else Vec2(0, -1) penetration oy contact_point Vec2( (a.position.x b.position.x) / 2.0, (a.position.y b.position.y) / 2.0 ) return normal, penetration, contact_point4.3 扩展到旋转矩形SAT 分离轴定理的前瞻如果哪一天你要支持矩形旋转Narrow Phase 里的矩形-矩形检测就必须升级为分离轴定理SAT。SAT 的核心思想是两个凸多边形不相交当且仅当存在一条轴使它们在这条轴上的投影不重叠可能的分离轴恰好是每个多边形边的法线方向。因此你只需要取两个矩形总共四条边的法线逐一做投影测试就能判断是否相交还能顺便求出最小穿透深度和法线。在 AABB 版本里检测轴正好就是 x 轴和 y 轴这就是为什么前面那段代码看起来这么简洁。将来升级 OBB 时只需把“绝对边界比较”换成“把顶点坐标旋转”再投影。核心逻辑完全复用。我建议你在备注里写下这一点方便以后翻代码时知道当初为什么要把 AABB 的先写出来。4.4 接触流形把几何信息打包传递给求解器每次 Narrow Phase 得到的信息我习惯用一个小结构体封装就叫Contactclass Contact: def __init__(self, body_a, body_b, normal, penetration, contact_point): self.a body_a self.b body_b self.normal normal self.penetration penetration self.point contact_point这里的关键点是normal的方向约定必须统一。我统一约定为从 A 指向 B这样冲量求解时只需要一个符号约定就能让两个物体分开。穿透深度是正值接触点用来计算切向摩擦和冲量产生转动的力臂。把所有这几项打包后后续的碰撞响应模块就可以完全不关心底层是圆还是矩形了——多态的效果。5. 碰撞响应冲量公式、恢复系数与摩擦力碰撞检测告诉我们“物体撞上了”碰撞响应则回答“接下来速度怎么变”。这里有一个基本物理概念上的分岔点为什么用“冲量”而不是“力”。因为碰撞发生的时间极短远小于一个模拟步长力在这个过程中不是一个有限可测的量真正起作用的是力对时间的积分也就是冲量。它直接改变的是速度而不是加速度。5.1 法向冲量的公式推导先假设两个物体只在接触点的法线方向上发生相互作用。定义法向相对速度vn dot(vB - vA, normal)。如果vn 0说明两个物体正在分离不需要响应如果vn 0说明它们在相互接近就必须施加冲量把速度改到分离方向。设冲量大小为j施加在 B 上的方向是normal施加在 A 上的方向是-normal。速度变化量为ΔvA -j * invMassA * normal ΔvB j * invMassB * normal碰撞后的法向相对速度应该满足恢复系数定义vn_after -e * vn_before把速度变化代入展开后解出j得到这个经典公式j -(1 e) * vn / (invMassA invMassB)注意分母是两者质量倒数之和这正是为什么第一章节坚持用invMass的原因——公式里根本不需要单独写质量。如果某个物体是静态的它的invMass 0分母自然只剩另一个物体的倒数效果上等价于撞了一面无限重的墙。落地代码如下。我把它写成一个独立的resolve_contact函数方便后续进行迭代求解时反复调用。def resolve_contact(contact, e): a, b contact.a, contact.b n contact.normal v_rel b.velocity - a.velocity vn v_rel.dot(n) if vn 0.0: return # 已经分离 inv_mass_sum a.inv_mass b.inv_mass if inv_mass_sum 1e-12: return # 两个物体都是静态 j -(1.0 e) * vn / inv_mass_sum impulse n * j a.velocity - impulse * a.inv_mass b.velocity impulse * b.inv_masse是恢复系数范围通常在 0 到 1 之间。e 0代表完全非弹性碰撞后法向相对速度为零两个物体“粘”在一起e 1代表完全弹性碰撞前后法向相对速度大小不变。真实世界里的物体都夹在中间比如混凝土约等于 0.0橡胶球约 0.8。如果你做的是休闲游戏恢复系数给得低一些反而更容易稳定高恢复系数会让能量积聚堆叠场景很快会抖起来。5.2 切向摩擦力库仑模型的简化实现法向冲量解决的是“弹开”切向冲量解决的是“滑动”。库仑摩擦模型里摩擦力的大小存在一个上限与法向力成正比这个比例系数就是摩擦系数 μ。摩擦系数分为静摩擦和动摩擦简易引擎里通常只用一个“动态摩擦系数”近似。计算切向相对速度的方法是把相对速度向量减去其在法线上的投影分量剩下的部分就是切向速度t v_rel - n * (v_rel.dot(n)) if t.length_sq() 1e-12: t t.normalized() else: return # 没有切向运动 vt v_rel.dot(t)切向冲量的目标是把切向相对速度尽可能降低但上限不能超过μ * |j|这里的j是前面求出的法向冲量绝对值。这样做的好处是摩擦不会反向推动物体产生不自然的运动。实现代码如下def apply_friction(contact, mu, j_normal): a, b contact.a, contact.b n contact.normal v_rel b.velocity - a.velocity t v_rel - n * v_rel.dot(n) if t.length_sq() 1e-12: return t t.normalized() vt v_rel.dot(t) inv_mass_sum a.inv_mass b.inv_mass jt -vt / inv_mass_sum max_friction mu * abs(j_normal) if jt max_friction: jt max_friction elif jt -max_friction: jt -max_friction impulse t * jt a.velocity - impulse * a.inv_mass b.velocity impulse * b.inv_massmu取值的物理参考冰面约 0.05橡胶与路面约 1.0。在 Demo 里我会用 0.3 左右既能看到滑行也不会滑得收不住。5.3 把法向冲量与摩擦合到一起的完整响应完整响应函数就是先算法向冲量再根据这个冲量计算摩擦。要注意顺序法向冲量改变了速度切向摩擦力基于的就是改变后的相对速度因此两者在一个接触点上是连续计算不是并行。把resolve_contact扩展成包含摩擦的版本后面迭代求解的每一轮都调用它def solve_contact(contact, e, mu): a, b contact.a, contact.b n contact.normal v_rel b.velocity - a.velocity vn v_rel.dot(n) if vn 0.0: return 0.0 inv_mass_sum a.inv_mass b.inv_mass j -(1.0 e) * vn / inv_mass_sum if inv_mass_sum 1e-12 else 0.0 impulse n * j a.velocity - impulse * a.inv_mass b.velocity impulse * b.inv_mass # 摩擦 v_rel2 b.velocity - a.velocity t v_rel2 - n * v_rel2.dot(n) if t.length_sq() 1e-12: t t.normalized() vt v_rel2.dot(t) jt -vt / inv_mass_sum if inv_mass_sum 1e-12 else 0.0 max_f mu * abs(j) if jt max_f: jt max_f elif jt -max_f: jt -max_f friction_impulse t * jt a.velocity - friction_impulse * a.inv_mass b.velocity friction_impulse * b.inv_mass return j到目前为止所有的冲量计算都还停留在“两个物体之间的一次性碰撞”阶段。一旦物体数量变多你会发现一个现象底下的方块被压得不断下沉上面的方块却还在乱弹。这就要进入下一章了。6. 稳定堆叠迭代求解与穿透修正我敢说每个写过简易物理引擎的人第一次把一堆方块叠起来时都经历过“快乐三秒然后全面崩塌”。原因是单个接触的冲量求解只顾局部没有协调全局。打个比方一个人站在一排人肩膀上最底下的人受到的力并不是最顶上那个人直接压下来的力而是整条压力链传导下来的总和只求解一次接触根本传递不到最底层。6.1 为什么单次求解会让堆叠崩溃假设三个方块 A、B、C 从上到下叠在一起C 在最底A 在最上。第一轮碰撞检测生成了两个接触对A-B 对和 B-C 对。如果只对每个接触对求解一次A 撞到 B 时 B 获得向下的速度增量但 B 和 C 之间的接触在检测时发现 B 的速度已经是穿透方向冲量会把它推回去一部分。这一套流程虽然能让每个局部碰撞看起来合理但整体上 A 的重力势能还没有“传导”到底部B 和 C 之间的接触力就不足以撑住 A于是 B 慢慢下压C 慢慢被压穿几秒后全塌。解决办法就是迭代。在同一帧内对所有接触对反复求解多次每次求解都会把上一轮产生的速度变化反馈到相邻接触相当于在压力链上多次传递信息。Box2D 这类引擎里那个velocityIterations参数就是干这个用的常见取值在 8 到 20 之间。迭代求解非常朴素在一个循环里反复调用solve_contact即可def solve_all_contacts(contacts, iterations10, e0.3, mu0.3): for _ in range(iterations): for contact in contacts: solve_contact(contact, e, mu)第一次迭代后问题已经减轻第二次之后压入深度会明显收敛五六次后肉眼基本分辨不出堆叠下沉。代价是计算量线性增加但对我们几十个物体的场景来说完全可以忽略。这个“简单到有点无脑”的方法恰恰是工业级引擎里约束求解器的雏形——只不过工业级引擎加入了矩阵求解和更复杂的松弛策略。6.2 穿透修正为什么不能只用冲量即使做了迭代冲量求解也只能修正速度已经发生的位移穿透并不会被完美拉回。想象一下一个球陷入了地面 0.1 像素冲量会让它的速度变为向上但下一秒它仍然会带着那个负的位移重新落入地面表现为肉眼可见的“地面之下抖动”。常见的修正手段叫投影法Projection直接把物体的位置沿法线方向推开穿透距离。但不能推满否则物体会像被弹簧弹走一样弹出老远。我用的策略是设定slop允许的微小穿透低于这个值不修正避免物体在稳定堆叠时一直微调抖动。slop通常取0.005。设定修正比例percent我习惯用0.8。一次修正 80% 的超出slop的穿透量留下来的微小穿透让自然碰撞去缓慢归位。def positional_correction(contacts, percent0.8, slop0.005): for contact in contacts: a, b contact.a, contact.b inv_mass_sum a.inv_mass b.inv_mass if inv_mass_sum 1e-12: continue correction_amount max(contact.penetration - slop, 0.0) / inv_mass_sum correction contact.normal * (correction_amount * percent) a.position - correction * a.inv_mass b.position correction * b.inv_mass注意修正的是位置而不是速度这本质上是在“作弊”——把物理上不该发生的穿透直接推回去。这个方法虽然不物理但在游戏物理中完全够用而且远比真正物理的“约束力求解”简单稳定。6.3 迭代次数与主动休眠迭代次数多到一定程度后效果增长非常有限反而可能引入新的抖动。经验上 8 到 12 次已经是休闲游戏的安全区间如果堆叠超过五层可以考虑扩容到 20 次。省不省这笔计算取决于你的目标平台——如果是移动端每一层循环都在烧 CPU值得用心调参。另一个节省计算的手段是让静止物体进入休眠。休眠思想很简单如果一个物体速度低于某个阈值比如 0.01且持续了若干帧就停止更新它只有当另一个活跃物体撞上它或者外力触发时才唤醒。休眠能显著减少堆叠场景的计算量因为静止堆放弃后不再需要反复求解接触。简易引擎里可以先不做但思路要留作扩展点。6.4 从“抖”到“稳”的调参顺序总结遇到堆叠不稳定不要瞎调恢复系数和摩擦系数。先按这个顺序检查物理步长是否固定不固定先修步长。迭代次数是否大于等于 8没有加上。穿透修正是否做了没有按percent0.8, slop0.005加上。恢复系数是否太高建议先全部设为 0.1 ~ 0.3。摩擦系数是否太低先试 0.5 再看效果。这套顺序我踩坑踩出来的照着调一遍大部分抖动问题能解决八成以上。7. 一个不加任何图形库的可运行 Demo现在把所有模块拼起来我写一个不依赖任何第三方库的完整脚本。它模拟了这样一幕几只球和几个方块从半空掉下来落在一片静止底板上最终堆叠成一堆。每帧输出模拟数据你可以直接跑也可以把它封装进自己的渲染循环里。这个脚本刻意不画画面因为画出图形会引入平台相关的依赖而物理引擎的核心逻辑是纯数值的。你完全可以用 HTML Canvas、Pygame、SDL 或者其它任何方式对接输出坐标数据。7.1 完整代码物理核心与模拟场景下面这段代码总计不到两百行就是一个可运行的最小物理引擎。我把前面所有讲解过的模块都整合进来了注释里标明了每一段对应文章里的哪一章。import random import math # ---------- Vec2 ---------- class Vec2: def __init__(self, x0.0, y0.0): self.x float(x) self.y float(y) def __add__(self, o): return Vec2(self.x o.x, self.y o.y) def __sub__(self, o): return Vec2(self.x - o.x, self.y - o.y) def __mul__(self, s): return Vec2(self.x * s, self.y * s) __rmul__ __mul__ def dot(self, o): return self.x * o.x self.y * o.y def cross(self, o): return self.x * o.y - self.y * o.x def length_sq(self): return self.x * self.x self.y * self.y def length(self): return math.sqrt(self.length_sq()) def normalized(self): l self.length() if l 1e-10: return Vec2(1, 0) inv 1.0 / l return Vec2(self.x * inv, self.y * inv) def perp(self): return Vec2(-self.y, self.x) # ---------- Body ---------- class Body: def __init__(self, shape, position, shape_data, mass1.0, restitution0.3): self.shape shape # circle / aabb self.shape_data shape_data # 半径或 Vec2(半宽, 半高) self.position position self.velocity Vec2() self.angle 0.0 self.angular_vel 0.0 self.force_acc Vec2() self.inv_mass 0.0 if mass 0.0 else 1.0 / mass self.inv_i 0.0 # 本Demo不使用旋转 self.restitution restitution def apply_force(self, f): self.force_acc f def step(self, gravity, dt): accel self.force_acc * self.inv_mass self.velocity accel * dt self.velocity gravity * dt self.position self.velocity * dt self.force_acc Vec2() # ---------- 碰撞检测 ---------- def circle_circle(body_a, body_b): delta body_b.position - body_a.position r_sum body_a.shape_data body_b.shape_data if delta.length_sq() r_sum * r_sum: return None dist delta.length() n delta * (1.0 / dist) if dist 1e-6 else Vec2(1, 0) pen r_sum - dist cp body_a.position n * (body_a.shape_data - pen * 0.5) return n, pen, cp def aabb_aabb(body_a, body_b): a_hw body_a.shape_data.x a_hh body_a.shape_data.y b_hw body_b.shape_data.x b_hh body_b.shape_data.y dx abs(body_b.position.x - body_a.position.x) dy abs(body_b.position.y - body_a.position.y) ox a_hw b_hw - dx oy a_hh b_hh - dy if ox 0.0 or oy 0.0: return None if ox oy: n Vec2(1, 0) if body_b.position.x body_a.position.x else Vec2(-1, 0) pen ox else: n Vec2(0, 1) if body_b.position.y body_a.position.y else Vec2(0, -1) pen oy cp Vec2( (body_a.position.x body_b.position.x) / 2.0, (body_a.position.y body_b.position.y) / 2.0 ) return n, pen, cp def circle_aabb(body_circle, body_rect): rect_center body_rect.position half_w body_rect.shape_data.x half_h body_rect.shape_data.y # 将圆心限制到矩形内部 closest Vec2( max(rect_center.x - half_w, min(body_circle.position.x, rect_center.x half_w)), max(rect_center.y - half_h, min(body_circle.position.y, rect_center.y half_h)) ) delta body_circle.position - closest if delta.length_sq() body_circle.shape_data * body_circle.shape_data: return None if delta.length_sq() 1e-12: # 圆心在矩形内部用最近的面推出 overlap_left body_circle.position.x - (rect_center.x - half_w) overlap_right (rect_center.x half_w) - body_circle.position.x overlap_bottom body_circle.position.y - (rect_center.y - half_h) overlap_top (rect_center.y half_h) - body_circle.position.y min_ov min(overlap_left, overlap_right, overlap_bottom, overlap_top) if min_ov overlap_left: n Vec2(-1, 0) elif min_ov overlap_right: n Vec2(1, 0) elif min_ov overlap_bottom: n Vec2(0, -1) else: n Vec2(0, 1) pen min_ov cp body_circle.position n * pen return n, pen, cp dist delta.length() n delta * (1.0 / dist) pen body_circle.shape_data - dist cp closest n * (body_circle.shape_data - pen * 0.5) return n, pen, cp def detect_contacts(bodies): contacts [] n len(bodies) for i in range(n): for j in range(i 1, n): a, b bodies[i], bodies[j] if a.inv_mass 0.0 and b.inv_mass 0.0: continue if a.shape circle and b.shape circle: result circle_circle(a, b) elif a.shape aabb and b.shape aabb: result aabb_aabb(a, b) else: if a.shape circle: result circle_aabb(a, b) else: # 交换顺序让 circle_aabb 始终接收 (圆, 矩形)法线取反 result circle_aabb(b, a) if result: n_vec, pen, cp result result (n_vec * -1.0, pen, cp) if result: normal, penetration, cp result # 注意检测阶段不改变任何速度 contacts.append((a, b, normal, penetration, cp)) return contacts # ---------- 碰撞响应 ---------- def solve_contact(a, b, normal, penetration, cp, e, mu): v_rel b.velocity - a.velocity vn v_rel.dot(normal) if vn 0.0: return 0.0 inv_sum a.inv_mass b.inv_mass if inv_sum 1e-12: return 0.0 j -(1.0 e) * vn / inv_sum impulse normal * j a.velocity - impulse * a.inv_mass b.velocity impulse * b.inv_mass # 摩擦 v_rel2 b.velocity - a.velocity tangent v_rel2 - normal * v_rel2.dot(normal) if tangent.length_sq() 1e-10: tangent tangent.normalized() vt v_rel2.dot(tangent) jt -vt / inv_sum max_friction mu * abs(j) if jt max_friction: jt max_friction elif jt -max_friction: jt -max_friction f_impulse tangent * jt a.velocity - f_impulse * a.inv_mass b.velocity f_impulse * b.inv_mass return j def solve_contacts(contacts, e, mu, iterations10): for _ in range(iterations): for contact in contacts: solve_contact(contact[0], contact[1], contact[2], contact[3], contact[4], e, mu) def positional_correction(contacts, percent0.8, slop0.005): for contact in contacts: a, b contact[0], contact[1] normal contact[2] pen contact[3] inv_sum a.inv_mass b.inv_mass if inv_sum 1e-12: continue amount max(pen - slop, 0.0) / inv_sum * percent corr normal * amount a.position - corr * a.inv_mass b.position corr * b.inv_mass # ---------- 主模拟 ---------- def physics_step(bodies, gravity, dt, e0.3, mu0.3, iterations10): # 1. 积分力与速度 for body in bodies: body.step(gravity, dt) # 2. 碰撞检测 contacts detect_contacts(bodies) # 3. 迭代求解冲量与摩擦 solve_contacts(contacts, e, mu, iterations) # 4. 位置修正 positional_correction(contacts) def build_scene(): bodies [] # 静态底板一个巨大的AABB块 ground Body(aabb, Vec2(0, -3.0), Vec2(10.0, 0.5), mass0.0) bodies.append(ground) # 动态方块 for i in range(4): bodies.append(Body(aabb, Vec2(-1.0 i * 0.7, 5.0 i * 0.2), Vec2(0.35, 0.35), mass1.0, restitution0.1)) # 动态圆球 for i in range(4): bodies.append(Body(circle, Vec2(2.0, 4.0 i * 0.5), 0.3, mass1.0, restitution0.3)) return bodies if __name__ __main__: random.seed(0) bodies build_scene() gravity Vec2(0.0, -9.81) dt 1.0 / 60.0 for frame in range(180): physics_step(bodies, gravity, dt, e0.2, mu0.4, iterations10) # 每30帧打印一次所有物体的位置 if frame % 30 0: info [] for b in bodies: if b.inv_mass 0.0: info.append(f({b.position.x:.2f}, {b.position.y:.2f})) print(fframe{frame:3d} objects: { .join(info)})把这段代码保存成.py文件直接运行你会看到 180 帧模拟 3 秒内所有动态物体从半空掉落、碰撞、最终趋于静止的过程。如果某一次运行出现了堆叠下沉检查是不是把迭代次数调得太低了正常跑下来最后输出的一行底部方块的位置应该稳定在y2.15上下静态地板中心在 -3.0半高 0.5顶面在 -2.5方块半高 0.35所以方块底面在 -2.5、中心应在 2.15。这个数值就是静置状态的标尺。7.2 怎样把数据变成可见画面控制台输出只能验证引擎数值但物理引擎的乐趣要看见才完整。如果你已经有一个渲染环境把bodies里每个Body的位置和形状画出来即可。比如用 HTML Canvas画一个圆就调用一次arc(x, y, r)画一个矩形就调用rect(x - half_w, y - half_h, half_w * 2, half_h * 2)。没有任何其它魔法物理核心永远只输出刚体状态数据。如果你偏好 Python那就在环境里安装 Pygame写一个几十行的绘制循环或者干脆把数据导出成 CSV用任何电子表格软件画散点图。核心逻辑和画面完全解耦这是这个设计最舒服的地方。8. 再往前一步从简易引擎到可商用引擎的路写到这里这套引擎已经能模拟基础刚体世界了但它离一个“能扛事的引擎”还有相当距离。我最后一次打开自己的这个项目时在代码旁边列了一张“还缺什么”的清单应该对你有同样的参考价值。8.1 下一步最值得做的四个扩展扩展一旋转刚体。给刚体加入角速度、角加速度并把接触点偏移量写进冲量公式。核心变化是冲量分母多出两项(rA × n)² * invIA (rB × n)² * invIB分别代表接触点冲量对两个物体转动的影响。扩展后你的引擎就支持了最基础的刚体动力学方块才能真正“翻滚”。扩展二SAT 支持任意凸多边形。轴对齐矩形版跑通后换成 SAT 前文已经讲过原理你可以慢慢把所有凸多边形统一到一个数据结构碰撞检测会更通用也能支持三角形。扩展三关节约束。距离关节、旋转关节、滑动关节本质上是把“接触约束”泛化成“更一般的约束”。理解冲量法后关节只是往约束矩阵里再塞几行约束方程。这类扩展是自制引擎最过瘾的部分因为你会突然拥有“设计机关”的能力。扩展四Broad Phase 换成空间哈希。当物体数量超过几十个时暴力遍历就开始拖后腿了。空间哈希每个格子里的物体维护一个链表每帧只检测当前格子及邻居格子里的候选对。它的实现难度不大但性能提升立竿见影。8.2 用个人经验做个收尾而不是替你做决定说句心里话现在常见的物理引擎都非常成熟如果你的目标是做一款完整商业游戏直接用成熟引擎是最快的路径。但我依然强烈建议每个想做物理相关功能的人亲手实现一遍简易引擎哪怕只是跑通这个 Demo。原因不是“造轮子锻炼能力”这种空话而是你在调堆叠参数、看它塌掉、再把它修好的过程中会建立起对物理引擎内部工作的直觉。这个直觉非常值钱它会在你以后用任何现成引擎时帮你快速定位问题绕开那些让新人崩溃的坑。如果你问我下一步怎么走我会说先把 Demo 加个简单的渲染层让它动起来感受一下自己亲手写出的“世界规则”。然后再给一个方块加上旋转你要过第一个坎是冲量分母那两项怎么加——加对后你会觉得整个世界都完整了。那一步的感受跟看到小球第一次不穿透地面时一样爽。

从零实现简易2D物理引擎：碰撞检测、冲量响应与稳定堆叠

相关推荐

YOLO训练过拟合怎么办？数据增强、正则化、早停三招解决

提示工程架构师如何驯服智能作曲：从情绪翻译到结构化指令

AI辅助开题报告写作：从框架到定稿的完整实操指南

更多资讯

Web版YOLO训练平台实战：FastAPI+SQLite单机部署全攻略

Agency Agents多智能体协作系统：四层架构解析与落地实践

从SEO到GEO：大模型时代的内容分发重构与模型引用优化实践

基于FastAPI和SQLite的YOLO训练网页化实践，单机部署方案详解

Python操作MySQL事务处理实战：从ACID到隔离级别与并发控制

SSMS使用指南：从连接实例到执行计划，避开掉五大常见坑

今日更新

电商智能运营体系从0到1：架构设计、落地步骤与避坑指南

压缩空气储能与热网热惯性联合优化调度：热电联产系统建模与求解

若依框架AI代码生成不稳定？用Agent Skill三层结构+五步法搞定

本周精选

金融AI落地指南：从调研报告到技术选型与工程实践

大模型上下文成本与注意力优化技术深度解析——为何O(n²)平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测

基于GMM的语音识别系统原理与MATLAB实现详解

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？