物体碰撞与摩擦的方法总结
本文禁止转载
B站：Heskey0
Contact and Friction Simulation for Computer Graphics(Siggraph course 2022)
相关的course：SIGGRAPH'20 Course: An Introduction to Physics-Based Animation
Syllabus
Basics：
- contact generation using discrete collision detection
- numerical techniques for solving the associated LCPs and NCPs
advanced topics:
- soft body contact approaches
- penalty and barrier functions
- anisotropic friction modeling
Chapter 1. Introduction to Contact Simulation
there are three main paradigms for contact simulation:
-
constraint based methods(很精确): constrained optimization（前三章）
- try to eliminate any penetration between bodies by applying an impulse at the point of contact
-
penalty based approaches: small abstract springs(也可以不使用弹簧) working between objects to keep them from overlapping. Hence, the springs "penalize" overlap.（第五章）
- similarities to penalty methods in numerical optimization: barrier methods
-
impulse-based methods: contact between objects is conceptualized as sequences of micro-collisions（本course notes不会提到，详情请见Bridson2002）
1.1 The Equation of Motion
Notations in this note:
- \(M\): mass matrix
- \(q\): position
- \(u\): velocity
- \(h=\Delta t\): timestep
1.2 Time Integration
Notations:
- \(u^-=u^n\): velocity of the last time step
- \(u^+=u^{n+1}\): velocity of the next time step
1.3 Constraints
two types of constraint equations:
- bilateral: \(\phi(q)=0\), hinges, ball-and-socket, and prismatic joints
- unilateral: \(\phi(q)\ge0\).
\(m\) constraint functions \(\phi(q)\in R^m\) implicitly define a manifold that is embedded in the \(n\)-dimensional space of the simulation degrees of freedom. We can instead formulate the constraint equations in terms of the velocities by computing the gradient of \(𝜙(q)\) with respect to \(q\), such that:
constraint equation:
- bilateral: \(Ju=0\).
- unilateral: \(Ju\ge0\).
约束就相当于沿 \(J\) 的方向施加一个力：the constraint forces of a system are computed as
where \(\hat \lambda\) are Lagrange multipliers
we need to integrate the constraint forces by applying the impulse \(hf_c\). Therefor, \(\lambda=h\hat\lambda\) is used to denote constraint impulse magnitudes
the equation of motion:
\(M\in R^{n\times n}\): the mass matrix
\(Mu\in R^n\): the momentum term(冲量)
\(f\): the applied force
we include the constraint impulse magnitudes as an implicit term, such that:
\(\lambda^+\in R^m\) is a vector of constraint impulse magnitudes
let us consider just the bilateral constraints in the simulation:
We use the first row to solve for \(u^+\) and substitute our result into the second row to obtain a reduced system that we can use to solve for \(\lambda^+\). This technique is called forming the Schur complement of the upper left block of the equation. It results in the reduced linear system:
i.e.
In the next section, we consider that contact can be modeled as a unilateral constraint with some special complementarity conditions on the constraint impulses and the relative velocities of bodies.
1.4 Non-interpenetration Contact Constraint
The state of this relationship between bodies is represented using a gap function, \(\phi\), which also happens to be the constraint function.
an impulse with magnitude \(\lambda_{\hat n}\) is applied at the contact location in order to keep the bodies from interpenetrating.
The direction of the impulse is determined by the constraint Jacobian \(J=[J_i\quad J_j]\).
当 \(\phi>0\) 时 \(\lambda_{\hat n}=0\)，当 \(\phi=0\) 时 \(\lambda_{\hat n}>0\)。这两个条件可以写成the position-level complementarity conditions of contact:
\(\perp\): the complementarity operator.
- if \(a,b\) are two scalars, then \(0\le a\quad\perp\quad b\ge0\) is defined to mean that \(a\ge0,b\ge0\), and \(ab=0\).
如果物体有比较大的互相远离的速度，则无需施加冲量。
the non-penetration complementary condition (互补条件):
\(v_{\hat n}=\dot\phi=\part\phi/\part t\): the normal component of the relative velocity between the bodies at a specific contact point.
1.4.1 Non-interpenetration Jacobian
the relative contact point velocity:
the relative velocity in the normal direction:
Using the shorthand \(r_A=(p-x_A)\) and \(r_B=(p-x_B)\), the definition of the skew-symmetric cross product matrix of a vector \(r=[x\quad y\quad z]^T\) is:
The relative velocity in the normal direction can then be written as the matrix-vector product:
\(v_A,v_B\) are the linear velocities of center of masses and \(\omega_A\) and \(\omega_B\) are the angular velocities. 如下图所示：
Left:
- the notation involved in computing the relative contact point velocity of two rigid bodies are depicted.
Right:
- The normal and vector arms and barycentric coordinates are the information needed to evaluate the contact Jacobian for two soft bodies modeled by tetrahedral elements.
1.4.2 Non-interpenetration Constraint for Soft Bodies
上一节的右图中
the point of contact \(p\):
where \(w_i,w_j,w_k,w_m\) are the barycentric coordinates of \(p\) with respect to the tetrahedron
1.5 The Coulomb Friction Law
there are two orthogonal unit vectors \(\hat t,\hat b\) that span the contact plane. the contact plane can be viewed as the shared tangent plane between two smooth surfaces at a point of contact \(p\). Orthogonal to the contact plane we have the contact normal direction \(\hat n\).
the contact space velocity:
\(\Delta v\): relative velocity of the surfaces in the world space in section 1.4.1
i.e.
Notations:
- \(\hat v_{\hat t}\): the tangential component of the contact space velocity (切平面上的速度)
- \(\lambda_{\hat t}\in R^2\): the tangential impulse (切平面上为阻止滑动产生的冲量)
1.5.1 Theory of Friction
Modern Robotics现代机器人学学习笔记12.2 - 摩擦锥
The exact isotropic planar Coulomb friction cone(摩擦锥) constraint can be written as: (收到的切向力不超出锥的范围，则物体不会移动)
\(\mu\): 摩擦系数，The Engineering ToolBox 中可以查阅到
The isotropic Coulomb friction law is in fact a two-part law:
- Slip: If the bodies are sliding relative to each other, then the direction of the friction force is opposed to the tangential relative velocity and its magnitude is \(\mu\lambda_{\hat n}\).
- Stick: If the objects are not moving relative to each other, then the friction force can have any direction. (前提是 \(||\lambda_{\hat t}||_2\le\mu\lambda_{\hat n}\))
We can express both cases mathematically as follows: (the stick-slip conditions)
第一式：摩擦力最大为 \(\mu F_N\), 其中 \(F_N\) 为正压力
第二式：滑动时，摩擦力达到最大值 \(\mu F_N\)
第三式：速度和摩擦力反向
This form of the isotropic Coulomb friction model is frequently used for modeling frictional contact, and it represents a non-linear complementarity problem(NCP) formulation of frictional contact
1.6 The Linear Complementarity Problem Model(LCP)
in this and the next section, we present details on how to linearize the Coulomb friction model into the standard LCP and BLCP forms. (linearize the 3D friction model using a polyhedral cone)
用棱锥近似圆锥，侧面面数设为 \(k\)：
用 \(k\) 个向量作为基，展开tangent space（positive span: 如下图：\(\hat t_1\)和\(\hat t_2\)相反，\(\hat t_3\)和\(\hat t_4\)相反，所以在表示空间中任意向量时，线性组合的系数均为非负数）
The positive span approximation allows us to rewrite \(||\lambda_{\hat t}||_2\le\mu\lambda_{\hat n}\) as:
where \(\lambda_{\hat t}=[\lambda_{\hat t_1}\quad...\quad\lambda_{\hat t_k}]^T\) is a vector of friction impulses along the tangent directions
If we let the maximum slipping velocity along any direction denoted by \(\beta\ge0\) then we must have:
where \(v_{\hat t}=[v_{\hat t_1}\quad...\quad v_{\hat t_k}]\) is the sliding velocity measured along each tangent vector and where \(e\) is a \(k\) dimensional vector of ones.
\(\beta\) is an estimate of the maximum sliding velocity along the directions of the positive span of tangent vectors.
Finally, we can write up the complete mathematical model for our derivation. Letting \(\lambda_{\hat t}=[\lambda_{\hat t_1}\quad...\quad\lambda_{\hat t_k}]\) be a vector of friction impulses, the linearized model, including complementarity conditions, can be formally stated as:
第一式：non-penetration constraint
第二式：in case we do have friction \(\lambda_{\hat t_i}>0\) for some \(i\), then \(\beta\) will estimate the maximum sliding velocity along the \(\hat t_i\)'s direction
第三式：the friction force is bounded by the Coulomb friction cone.
Jacobian: (For a general system consisting of \(p\) contacts)
Finally, the constrained equations of motion for the system can be written as:
where:
- \(\bar\mu=diag([\mu_1\quad...\quad\mu_p])\) is a diagonal matrix containing friction coefficients of each contact
- \(E=diag([e^T_1\quad...\quad e^T_p])\) is a block diagonal matrix containing \(k\) dimensional vectors of ones.
- \(s_{\hat n},s_\beta,s_\lambda\) helps to write the complementarity problem in a more concise way.
By applying the Schur complement “trick”, we can compute the reduced system:
\[(C+GM^{-1}G^T)\pmatrix{\lambda^+_{\hat n}\\ \lambda^+_{\hat t}\\ \beta }+GM^{-1}(Mu+hf)=\pmatrix{s_{\hat n}\\s_\beta\\s_\lambda} \]i.e.
\[Ax+b=\pmatrix{s_{\hat n}\\s_\beta\\s_\lambda} \]where:
\[C=\pmatrix{0&0&0\\0&0&E\\\bar\mu&-E^T&0},\textrm{and},G^T=\pmatrix{J^T_{\hat t}&J^T_{\hat n}&0} \]
The first equation can be written concisely and compactly in the form of a standard LCP:
This illustrates how powerful the LCP model really is. By simply assembling the \(A\) matrix and b vector, we can call our favorite LCP solver to compute \(x\).
1.7 The Boxed Linear Complementarity Problem Model(BLCP)
the friction cone approximated by a box, where the limits of friction impulses in each tangential direction \(t_i\) are determine by the box limit equation \(\lambda_{t_i}^{lo}\le\lambda_{t_i}\le\lambda_{t_i}^{hi}\)
discretizing the friction cone using two orthogonal tangent directions:
where:
- \(h_{\hat t_1}^{lo}=-\mu\lambda_{\hat n}\) and \(h_{\hat t_1}^{hi}=\mu\lambda_{\hat n}\)
Some iterative solvers update these friction bounds during the solve based on the magnitude of the normal forces, yielding a four-sided pyramidal cone that is larger than the true friction cone.
Problems in this more general class are called boxed linear complementarity problems (BLCP).
We note that by decomposing the residual velocity as \(v_i=^+v_i-^-v_i\), the BLCP can be written as three separate LCPs:
Jacobian:
the global Jacobian \(J=[J_1^T\quad...\quad J^T_p]^T\).
The BLCP of the constrained equations of motion:
i.e. \(Ax+b=[0\quad v]^T\)
Using Schur complement results in:（利用方程组的第一行解出 \(u^+\)，带入方程组第二行）
i.e. \(A\lambda^++b=v\)
1.8 The Cone Complementarity Problem Model(CCP)
1.9 Constraint Stabilization
1.10 Soft vs. Rigit Body
Chapter 2. Contact Generation
contact position is a location that represents the overlapping volume of the two bodies
contact normal is the direction of a restorative force that is applied to separate the two bodies
ideal touching state: 刚好接触不渗透
描述碰撞接触的三个重要属性：
-
Position: Positions are used as the point of action when applying contact forces. ghost torques may be introduced if \(p_A\ne p_B\).
- 接触点：\(p\)。无穿透时 \(p=p_A=p_B\)。有穿透时 \(p=\frac12(p_A+p_B)\)。
-
Normal: normal is computed and used to apply non-interpenetration forces in the correct direction.
- 通常 \(\hat n_A=-\hat n_B\). Use \(\hat n=\hat n_A\) or \(\hat n=\hat n_B\) as the normal associated with a contact point.
-
Penetration Measure: The measure is often used to add stabilization terms to counter drift errors.
- gap function: \(\phi\). the distance one would need to move along \(\hat n\) to bring the two objects into a touching state
2.1 Analytic Shapes
可以先对复杂物体的简单包围体进行碰撞检测（sphere, capsules, box）
- Sphere: center point \(c\in R^3\) and radius \(r\in R\).
- Box(OBB oriented bounding boxes): center point \(c\in R^3\), rotation matrix \(R\in SO(3)\) and half-width extents along the local axes \(h\in R^3\) (边长的一半)
2.1.1 Sphere-Sphere Intersection
if gap function is negative or zero, a contact is generated:
the unit length normal:
the contact point(the centor of the overlapping region):
2.1.2 Sphere-OBB Intersection
first transforming the center of the sphere into the local coordinate system of the box:
each coordinate \(i\in\{x,y,z\}\) of the closest point(contact point) \(g\) is computed as:
if \(||g-c_{\textrm{local}}||<r\), then there is an intersection.
special case(deep penetration): sphere center is entirely inside the box, i.e. \(g=c_{\textrm{local}}\).
(1) Shallow penetration
Shallow penetration: the transformed sphere center lies outside the box extents in at least one dimension.
The penetration \(\phi\) is computed as:
The normal: the unit vector that points from the sphere center to the closest point transformed in the global reference frame
(2) Deep penetration
the normal of the closest box face also determines the contact normal.
transform the closest point to the global coordinate system:
2.1.3 OBB-OBB Intersection
separating axis theorem (SAT): 将一束平行光照射到物体上，如果在任何角度的照射下，两个物体的影子都会重合，那么这两个物体一定相交（只需要在两个多边形所在的每个边上投影即可）
- \(d_{AB}=c_B-c_A\)，将两个OBB投影到轴 \(\hat l\).
- 由于由于OBB的对边是平行的，所以只需要三条投影轴（分别是OBB三个面的法向量）
- separation equation: \(s=|d_{AB}\cdot\hat l|-(|h_A\cdot\hat l|+|h_B\cdot\hat l|)\). If \(s\le0\) for all SAT tests, then an overlap exists and a contact will be generated.
2.2 Mesh-based Representation
BVH / spatial hashing: a broad phase collision detection method determines triangles that are sufficiently close and which pair-wise triangles need further processing
consider no penetration or separation between the two objects; they are simply touching: Corner points of polygonal contact areas can be represented by pairs for mesh features. The vertex-face \((V,F)\) and edge-edge \((E,E)\) combinations forms a sufficient set of pair-types.
2.2.1 CCD of Mesh Features
There are two flavors of continuous collision detection:
- look forward in time, from the current instance of invocation to the next expected instance
- look backward in time, from the current instance of invocation to the previous instance.
Consider three vertices of a triangle move with velocities \(u_0,u_1,u_2\), note as \(u_i\).
- for a forward prediction scheme, the candidate position: \(x^\prime_i=x_i+hu_i\).
- for a backward prediction scheme, the candidate position: \(x^\prime_i=x_i-hu_i\).
检测步骤：
-
A simple preliminary test would be to create two axis aligned bounding boxes (AABBs) that enclose the triangles– one using the current positions and the other using the candidate positions– and test them for overlap. Only if the bounding boxes overlap can there be a contact that occurred between the mesh features during the interval \(h\).
-
The next step is to examine all possible contact types for the triangle pair.
2.2.2 CCD of Face-Vertex Pair
face \(x_1,x_2,x_3\), vertex \(x_4\), 发生接触的方程为 (Co-planar test):
This is a cubic polynomial in \(t\).
碰撞接触的三个重要属性：
- contact point: \(p=x_4+u_4t\)
- normal: \(n=(x_2(t)-x_1(t))\times(x_3(t)-x_1(t))\)
- penetration depth: \(\phi=0\)
2.2.3 CCD of Edge-Edge Pair
edge \(x_1,x_2\), edge \(x_3,x_4\),
- firse, test whether the edges are parallel:
- compute the closest points between the two lines using normal equation(Bridson 2002中有提到)
碰撞接触的三个重要属性：
- contact point: \(p=(x_1-x_3)+a(x_2-x_1)-b(x_4-x_3)\)
- normal: edge cross product \(\hat n=e_A\times e_B/||e_A\times e_B||\)
- penetration depth: \(\phi=0\)
2.2.4 Multivariate Style of CCD
(TOG 2021) A Large-Scale Benchmark and an Inclusion-Based Algorithm for Continuous Collision Detection
For V-F case:
\(u,v\) are the barycentric coordinates of the triangle
\(t\in[0,h]\) and \(u\ge0,v\ge0\).
For E-E case:
\(t\in[0,h]\) and \(a,b\in[0,1]\).
2.3 SDF
SDF: \(s(x)\), 内负外正
2.3.1 SDF-Point Intersection
Test if \(s(x)\le0\) is true
gap function: \(\phi=s(x)\)
contact normal: \(\hat n=\nabla s(x)^T\)
contact point: two choices:
- \(p=x-\nabla s(x)^Ts(x)\)
- \(p=x\)
2.3.2 SDF-Mesh Intersection
(Macklin 2020) Local Optimization for Robust Signed Distance Field Collision
将Mesh B的顶点转换到Mesh A的local coordinate：
where \(c_A ∈ R^3 , R_A ∈ 𝑆𝑂(3)\) and \(c_B ∈ R^3 , R_B ∈ 𝑆𝑂(3)\) are the positions and rotation matrices of bodies A and B respectively
逐个顶点做SDF-Point Intersection test，得到\(\phi,\hat n,p\)
Chapter 3. Numerical Methods
In this section, we present a collection of numerical methods for solving the frictional contact problems.
两类经典方法：
- pivoting methods: focus on determining a labeling of variables based on the complementarity conditions. Once the labeling is known, an exact solution to the complementarity problem may be computed if a direct solver is used to solve the resulting subproblem.
- iterative methods: make incremental progress toward the exact solution, with better approximations being obtained at each iteration. Additionally, there are hybrid approaches that combine these schemes.
This table presents an overview of the main numerical methods covered in this chapter as well as a summary of some their traits.
3.1 Pivoting Methods
Pivoting methods exploit the combinatorial nature of complementarity problems to find a labeling of the variables(define index set) based on the complementarity and feasibility conditions imposed by the model. With this labeling, which is sometimes referred to as the index set, it is possible to identify constraint variables as being free or tight, and the latter case is further decomposed into tight lower and tight upper variables.
Notations:
- \(F\): the set of free variables
- \(L\): the set of tight lower variables
- \(U\): the set of tight upper variables
- \(A_{F,L}\): sub-block of \(A\) contains row indices from \(F\) and column indices from \(L\).
variables in \(i\in F\) are those for which the value of \(x_i\) is not restricted by the feasibility condition of the constraint. Whereas tight variables \(i\in L\cup U\) have a value which is determined by the feasibility conditions
To better understand how this partitioning works, let us consider the BLCP formulation and define the index sets:
举个栗子：
This labeling permits a partitioning of a linear system based on the index set, For instance:
对于这个方程组的第一行，由于 \(\forall i\in F:v_i=0\)，第一行就可以写为：
Hence, computing \(x_F\) simply requires solving an unconstrained linear system.
The objective of pivoting methods is therefore to correctly determine the index set of constraint variables that solves the complementarity problem.
3.1.1 Incremental Pivoting
(SIG 1994)(Baraff) Fast contact force computation for nonpenetrating rigid bodies.
3.1.2 Principal Pivoting Methods
3.1.3 Block Principal Pivoting
Chapter 4. Soft Body Contact
Chapter 5. Penalty and Barrier Methods
本章介绍penalty-based methods和barrier methods
penalty-based and constraint-based approaches are frequently combined for many implementations, giving a sort of “hybrid” simulator
5.1 Overview of Classical Penalty Methods
a spring-damper system is used for penalizing penetrations
The central idea behind penalty methods is to find a force \(f\) and torque \(\tau\) at a given instant that can be included in the dynamical equations of motion in order to perform the numerical integration.
In a practical implementation, one would iterate over contact points and compute the penalty force only once per contact and then distribute the spring force to the two bodies that are in contact. In its most pure form the simulation loop of the penalty method can now be summarized as:
- Detect contact points (run collision detection)
- Compute and accumulate spring forces
- Integrate equations of motion forward in time
the penalty force:
where:
- the \(k\)th contact point has the contact normal \(n_k\).
- penetration depth \(\phi_k\).
- \(v_k\): relative contact velocity at \(k\)th point
- \(k\): stiffness
- \(b\): damping
5.2 Barrier Methods
Penalty-based methods are limited by the choice of the penalty stiffness \(k\). Beyond this, no matter the choice of \(k\) there always exists a large enough velocity between contacting bodies that will result in tunneling. To address this problem we can turn to barrier methods.
Barrier methods are similar to penalty-based methods but apply penalties that diverge as objects get closer. (This is done through the use of so called barrier functions, which grows to infinity as our gap shrinks to zero.)
5.2.1 Asynchronous Contact Mechanics
a simple quadratic penalty:
where:
- \(\eta\): surface thickness
- \(x_a,x_b\): contact point
use this penalty to define a contact potential
where:
- \(k\): the penalty stiffness
5.3 Incremental Potential Contact
\(d\) 通常用unsigned distance表示
5.3.1 Consistent Distance Function
distance between point-triangle:
distance between edge-edge:
both of these constrained minimizations have a closed form solution.
这一部分的代码：ipc-toolkit/src/distance at main · ipc-sim/ipc-toolkit (github.com)
5.3.2 Optimization-based Time Integration
objective function:
可以使用gradient descent或者Newton's method做优化
蓝色的 PDProject 表示将Hessian 投影到 Positive Definite（这个操作是为了让newton不去找最近的鞍点，而去找极小点）
While Newton does not promise global convergence, the benefit of time dependent simulation is the previous time step usually provides a good warm start.
正定化：
Overview: The projection to PD is done by first projecting the elasticity, barrier, and friction Hessians to positive semi-definite (PSD). Then when we add the lumped mass matrix (Hessian of the inertia term in \(f\) ) which is a PD diagonal matrix, we get a resulting matrix \(H\) that is PD.
Projecting a matrix to PSD can be done by computing its Eigendecomposition and then setting all negative Eigenvalues to zero. However, computing a full Eigendecomposition can be prohibitively expensive. To avoid this large cost, we can instead project the local Hessian blocks rather than the entire matrix. This is, when we compute the Hessian of the \(k\) th contact term (which is a small \(12\times12\) matrix at most) we project it to PSD and then assemble the full matrix as usual. The same can be done for the elasticity Hessian on a per element bases.
Limitation: projecting to PD guarantees a descent direction it can worsen the rate of convergence of the method.
5.3.3 CCD-Aware Filtered Line-Search
A critical part of the IPC algorithm is ensuring the optimization does not take a step that jumps over the barrier. To accomplish this CCD is added to the line-search within the optimization.
line-search: finding a step length \(\alpha\in[0,,1]\) that decreases the energy between \(q_i\) and \(q_i+\Delta q_i\).
backtracking line search: \(\alpha\) is initialized to one. Then \(\alpha\) is repeatedly halved until a descent step \((f(q_i+\alpha\Delta q_i)<f(q_i))\) is found.
CCD-aware backtracking line-search: the initial value of \(\alpha\) is set to the minimum time of impact. By doing this it is guaranteed that all \(\alpha\le\textrm{CCD}(q,q+\Delta x)\) are intersection free.
5.3.4 Variational Friction Model
we can defines friction forces variationally by maximizing the dissipation rate subject to the Coulomb constraint.
challenge when integrating such forces: the transition between sticking and sliding comes with a discontinuous jump in both the force magnitude and direction. To address this, the transition between static and dynamic friction is mollified using a function:
where:
- \(s\): is the input (speed)
- \(\epsilon_v\): controls the accuracy of the model with smaller values approaching the full discontinuous transition(similar to \(\hat d\) in the contact formulation)
This lead to a smooth friction force:
where:
- \(T_k(q)\in R^{3n\times2}\): sliding basis at the point of contact
- \(N_k(q)\): the magnitude of the normal force
- \(\mu\): local friction coefficient
- \(\Delta\tilde q_k=T_k(q)(q-q^t)\in R^2\): the relative sliding displacement of the \(k\)th contacting point over the timestep.
By setting \(N_k\) and \(T_k\) to some fixed value \(N_k^{\textrm{lagged}},T_k^{\textrm{lagged}}\) we can define a dissipative potential as
where:
- \(M^\prime_{\epsilon_v}(s)=m_{\epsilon_v}(s)\)
the friction force is related to this potential by \(F_k(x)=-\nabla D_k(x)\).
This dissipative potential can then be added to the objective function. Again this leads to an unconstrained optimization that can be solved with Newton’s method using a CCD-aware line-search.
