游戲或者仿真中要提供接近於真實世界的完整觸覺反饋體驗，需要VR頭戴設備、控制器、外骨骼甚至是行走模擬裝置的配合。然而，人類的觸覺系統極其敏感，普通人打麻將就能用手指輕輕松松地摸出牌面。在目前的技術基礎上，機器很難還原真實的人類觸覺反饋，能做的只是在特定的內容和場景中盡量滿足用戶的反饋體驗。比如在VR游戲中抓取武器和物品時，幾厘米的偏差，沒有准確還原物體材質和紋路，對用戶的實際體驗影響並不大。除了高精度的外設，要實現觸覺或力反饋，還需要強大的物理引擎，虛擬世界中的復雜物體的建模、與用戶肢體的觸碰非常消耗計算資源。
- 牛頓運動定律與物理引擎
- 慣性定律：如果物體沒有受到外力作用，則運動狀態（靜止或勻速直線運動）不會發生改變；
- 牛頓第二定律：物體的加速度跟其所受合外力成正比，跟物體質量成反比，加速度方向與合外力方向相同，即$F=ma$；
- 牛頓第三定律：力是物體與物體之間的相互作用，作用力與反作用力總是大小相等方向相反，作用在同一直線上。
根據牛頓三大運動定律，我們可以創造物理引擎來使虛擬場景中的物體產生動態行為，提供沉浸式體驗。對於一個交互式游戲來說其基本運行流程如下：
計算物體間的碰撞是物理引擎最核心的部分。一般來說，先要解決如何高效地檢測到碰撞的產生（碰撞檢測），以及如何確定碰撞點及方向，之后我們就可以求得碰撞體的受力情況，從而根據牛頓運動定律計算出它們將要產生的平動和轉動。最后將場景中物體的位置和姿態輸出給圖形引擎去渲染。
但如何求得力和力矩？這便是一個很復雜的問題了。像重力可以直接影響物體所受力，摩擦可以直接影響物體所受力矩，這些都很簡單。較為復雜的就是兩個物體間的碰撞了，物體引擎中有一半以上的代碼是用來計算碰撞的。
在兩物體接觸面上存在着接觸面法向力（沖擊力或接觸力）和接觸面切向力（摩擦力）兩種作用。當接觸不連續時產生的接觸法向力就是沖擊力。由於摩擦力本身也是很復雜的問題，很多仿真軟件中對它作了簡化，采用比較簡單的庫倫摩擦力模型來計算：$$F_f=\mu F_n$$
- 用於碰撞的動力學基本定理
由於碰撞時間極短，通常只有千分之一甚至萬分之一秒，因此所產生的力非常巨大。這種產生在碰撞中，作用時間極短，數值巨大的力稱為碰撞力瞬時力。瞬時力的沖量稱之為碰撞沖量。瞬時力不僅數值巨大，而且隨時間迅速變化，其規律非常復雜，難以確定。碰撞過程中除了由碰撞力引起物體塑性變形外，同時還伴隨着發聲、發光和發熱等機械能轉換為其它形式能量的現象。因此，在研究碰撞問題時，一般並不去討論瞬時力本身，而只討論它的沖量及產生的總效果。研究碰撞問題，各微分形式的動力學基本定理不能直接應用，一般用積分形式的動量定理和動量矩定理。
當兩物體碰撞時，過接觸點作物體表面的公法線，此法線稱為碰撞法線。若兩物體的質心均位於碰撞接觸點的公法線上，則為對心碰撞；否則為偏心碰撞。
1. 沖量定理
對於由$n$個質點組成的質點系，其中第$i$個質點受到的碰撞沖量可分為外碰撞沖量$\bf I_i^{(e)}$和內碰撞沖量$\bf I_i^{(i)}$，由於內碰撞沖量大小相等方向相反，則質點系在碰撞開始和結束時的動量的改變等於作用於質點系的外碰撞沖量的矢量和：$$ m{\bf u_c} - m \bf v_c=\sum I_i^{(e)}$$
其中，$m$為質點系質量，$\bf v_c$和$\bf u_c$分別為質心碰撞開始和結束時的速度。
2. 沖量矩定理
質點系對質心的動量矩的改變等於外碰撞沖量對質心矩的矢量和：$$ \bf L_{c2}-L_{c1}=\sum M_c (I_i^{(e)})$$
式中，$\bf L_{c1}$和$\bf L_{c2}$分別表示碰撞開始和結束時質點系對質心的動量矩矢；$\sum M_c (I_i^{(e)})$為外碰撞沖量對質心之矩的矢量和。
對於平面運動剛體的碰撞問題，可應用質點系相對於質心的沖量矩定理來描述轉動部分。根據剛體平面運動特征，相對於質心的沖量矩定理在其平面內視為代數量，即有$$L_c=J_c \omega$$
式子中$J_c$為剛體對過質心且與其對稱平面垂直的軸的轉動慣量，$\omega$為剛體轉動的角度。由沖量矩定理有：$$L_{c2}-L_{c1}=J_c \omega_2-J_c \omega_1=\sum M_c (I_i^{(e)})$$
式中，$\omega_1$和$\omega_2$分別為平面運動剛體碰撞前、后的角速度。
- 恢復系數
碰撞過程可分為兩個階段。開始碰撞到物體速度為零的過程為第一階段，稱為變形階段。碰撞沖量為$\bf I_1$，則根據沖量定理：$$0-(-mv)=\it I_1$$
之后，物體恢復彈性變形到碰撞結束的過程為第二階段，稱為變形恢復階段。設碰撞沖量為$\bf I_2$，則根據沖量定理：$$mu-0=\it I_2$$
經研究發現，對於材料確定的物體發生正碰撞時，碰撞前后的速度大小之比幾乎是不變的，等於一個常數$e$，即$$e=\frac{u}{v}=\frac{I_2}{I_1}$$
$e$稱為恢復系數。恢復系數$e$表示物體在碰撞前后速度的恢復程度和物體變形的恢復程度，也反映了碰撞中機械能損失的程度。$e$越小，動能損失越大；反之，動能損失越小。一般情況下，$u<v$，恢復系數小於1。但理論上，恢復系數也可以大於1。例如當兩個手雷碰撞在一起產生爆炸，化學能轉換為機械能（能量增加），碰撞后的速度大於碰撞前速度。另外需要注意：恢復系數是兩個碰撞物體之間的共同性質，但在許多文獻中恢復系數常寫為單個物體固有屬性（沒有提這物體到底是與哪個物體相互碰撞），在這種情況下，第二個物體被假定為完全彈性剛體。
- ADAMS中的沖擊力模型
We can either use information about the depth of the collision and generate a very large force acting on the two objects (like a very stiff spring due to the “compression” of the object’s material as one penetrates the other), or we can use an impulse response which means we simply modify the objects’ momentum without any regard for the forces involved.
ADAMS中計算法向接觸力大小有兩種模型：沖擊函數模型（Impact）和恢復系數模型（Restitution）。恢復系數模型基於沖量理論。沖擊函數模型用一個彈簧-阻尼模型來表示（碰撞力由兩個部分組成：一個是由於兩個構件之間的相互切入而產生的彈性力；另一個是由於相對速度產生的阻尼力），沖擊函數表示為：$$F_n=k\cdot g^e+step(g,0,0,d_{max},c_{max}) \cdot\frac{dg}{dt}$$
式中：$g$——接觸物體之間的穿透深度；$F_n$——法向接觸力大小；$k$——剛度系數；$e$——碰撞系數，反映了材料的非線性程度；$c_{max}$——最大阻尼系數，表示物體碰撞時的能量損失；$d_{max}$——最大切入深度，它決定了何時阻尼力達到最大；為了防止碰撞過程中阻尼力的不連續，式中采用了step函數，其形式為$step(x,x_0,h_0,x_1,h_1)$
[Damping Coefficient versus Penetration]
- VREP中獲取接觸力的方法
機器人仿真中常常需要測量腳底壓力分布，或是手爪夾持力等接觸力。在VREP中獲取接觸力有兩種方法：一種是通過在末端添加力傳感器進行測量，另一種是直接利用物理引擎計算結果獲取接觸力信息，相關函數為simGetContactInfo：
使用simGetContactInfo函數可以直接獲取接觸力信息，不用添加力傳感器等物體，但是用法稍微復雜。VREP中默認情況下物理引擎運行速度是仿真速度的10倍，就是說如果仿真一步需要50ms，那么物理引擎已經運算了10步（5ms一步），如果dynamicPass參數設為sim_handle_all，則一步仿真將會返回10次接觸力的結果。如果要利用接觸力的信息，需要解析這些力屬於哪兩個接觸體。參考官方論壇中的帖子Get Friction between leg-tip and ground：
Is it possible to detect the friction between the tip of a robot-leg and the ground?
You have several ways of doing this:
- you can attach a force sensor at the tip of your leg, then attach a contact sphere to the force sensor (i.e. legTip --> forceSensor --> contactSphere). Then you can read the force/torque in the force/torque sensor
- you can parse all contact points and vectors for one simulation step and try to figure out which one(s) correspond to a specific leg/ground pair. This is more complicated, but you do not require an additional contact sphere as in above's first point. For that, have a look at the demo model Models/other/contact display.ttm (which uses the API function simGetContactInfo)
VREP中提供了一個contact display模型（Model browser --> other文件夾下），利用simGetContactInfo函數方便的為用戶顯示接觸點和接觸力的大小以及方向。修改其lua代碼，在狀態欄中顯示每次獲取的接觸力信息。從圖中可以看出立方體與地面接觸點有4個，將這些接觸點對應的接觸力加起來就是立方體重量。
核心代碼是下面這句，返回值objectsInContact是相互接觸兩物體的句柄，contactPt是接觸點的三維坐標，forceDirectionAndAmplitude是力向量的三維坐標，幅值代表大小：
objectsInContact,contactPt,forceDirectionAndAmplitude = simGetContactInfo(sim_handle_all, sim_handle_all, index)
仿真時可以看到四個點的接觸力數值是一直在跳動的（理論上如果接觸面是平面，四個點接觸力應該一樣），也許是三點確定一個平面，多一個冗余的接觸點就會導致計算不穩定？下圖中的球與地面只有一個接觸點，仿真時的接觸力要穩定許多。而且接觸力的仿真結果還與具體的物理引擎有關，Bullet和ODE的結果就比Vortex和Newton差很多...
為了得到穩定可靠地接觸力，可以對一次仿真中的所有接觸力求平均值（參考官方論壇帖子Impact force on collision），下面修改了contact display中的代碼，將每一步獲取的力（對於單個立方體與地面接觸模型，每一步仿真共有40個接觸點的接觸力）求均值，計算立方體重量：
if (sim_call_type==sim_childscriptcall_initialization) then black={0,0,0} purple={1,0,1} lightBlue={0,1,1} options=0 forceVectorScaling=simGetScriptSimulationParameter(sim_handle_self,'forceVectorScaling') forceVectorWidth=simGetScriptSimulationParameter(sim_handle_self,'forceVectorWidth') contactPointSize=simGetScriptSimulationParameter(sim_handle_self,'contactPointSize') overlayDisplay=simGetScriptSimulationParameter(sim_handle_self,'overlayDisplay') --sim_drawing_overlay: if specified, then items are drawn on top of other objects and are (almost) always visible if (overlayDisplay) then options=options+sim_drawing_overlay end -- Add a line and a sphere container: --sim_drawing_lines:items are pixel-sized lines. 6 values per item (x0,y0,z0,x1,y1,z1) lineContainer = simAddDrawingObject(sim_drawing_lines+options,forceVectorWidth,0,-1,1000,black,black,black,purple) --items are "sphere points". 3 values per item (x,y,z) sphereContainer = simAddDrawingObject(sim_drawing_spherepoints+options,contactPointSize,0,-1,1000,black,black,black,lightBlue) ObjectHandle = simGetObjectHandle("Cuboid") end if (sim_call_type==sim_childscriptcall_cleanup) then -- Remove the containers: simRemoveDrawingObject(lineContainer) simRemoveDrawingObject(sphereContainer) end if (sim_call_type==sim_childscriptcall_sensing) then -- empty the containers: simAddDrawingObjectItem(lineContainer,nil) simAddDrawingObjectItem(sphereContainer,nil) -- Fill the containers with contact information: index = 0 -- zero-based index of the contact to retrieve sum = 0 while (true) do -- Retrieves contact point information of a dynamic simulation pass objectsInContact,contactPt,forceDirectionAndAmplitude = simGetContactInfo(sim_handle_all, ObjectHandle, index) if (objectsInContact) then sum = sum + forceDirectionAndAmplitude[3] if (index == 39) then str = string.format("%f N", sum/10) simAddStatusbarMessage(str) -- Adds a message to the status bar end line={contactPt[1],contactPt[2],contactPt[3], 0,0,0} line[4]=contactPt[1]+forceDirectionAndAmplitude[1]*forceVectorScaling line[5]=contactPt[2]+forceDirectionAndAmplitude[2]*forceVectorScaling line[6]=contactPt[3]+forceDirectionAndAmplitude[3]*forceVectorScaling simAddDrawingObjectItem(lineContainer,line) -- Draw force vector simAddDrawingObjectItem(sphereContainer,{contactPt[1],contactPt[2],contactPt[3]}) -- Draw contact point index = index + 1 else break end end end
代碼中還涉及到了腳本仿真參數，通過設置仿真參數可以方便地修改控制腳本中的某些變量，比如線的寬度，接觸點大小，比例因子等等。The main script and each child script have a list of simulation parameters. Those parameters can be used as a quick way of adjusting values of a specific model. 在腳本文件中可以通過simGetScriptSimulationParameter函數來獲取用戶輸入的參數。仿真參數圖標位於腳本圖標右側，如下圖所示分別是沒有定義和定義了仿真參數的圖標樣式：
[Script simulation parameter icons (1) empty parameter list, (2) non-empty parameter list]
雙擊仿真參數圖標可以打開參數設置對話框，點擊Add new parameter按鈕可以插入新的仿真參數。Value處輸入參數的值，Unit為參數單位（這里只起到提示作用）
[Script simulation parameters dialog]
- Parameter is private: if enabled, then the selected parameter is not shown during a simulation (in that case the parameter is probably not meant to be modified during a simulation).
- Parameter is persistent: if enabled, then the selected parameter will not be restored to its original value at simulation end.
獲取仿真參數的函數simGetScriptSimulationParameter原型如下：
boolean/number/string parameterValue=simGetScriptSimulationParameter(number scriptHandle,string parameterName,boolean forceStringReturn=false)
函數參數scriptHandle為腳本句柄，或取值sim_handle_main_script或sim_handle_self；參數parameterName為仿真參數名（取值為腳本參數對話框中插入的參數名）；參數forceStringReturn默認為false，該參數用於指定是否將仿真參數值當作字符串。比如，如果該參數設為true，則contactPointSize的值0.01會被函數返回為string類型，即"0.01"，而非number類型。
下面是一個斜面接觸的例子，如下圖所示質量為1Kg的立方體由於靜摩擦力在30°斜面上保持靜止，可以看出下方兩個接觸點的接觸力要明顯大於上方的兩個（不像與水平面接觸時那樣均勻）。計算四個點接觸力矢量的幅值並相加可以求得立方體與斜面接觸力的大小，理論上這個力為mgcos30°≈8.66N，實際仿真時輸出的結果與理論值一致：
if (sim_call_type==sim_childscriptcall_initialization) then black={0,0,0} purple={1,0,1} lightBlue={0,1,1} options=0 forceVectorScaling=simGetScriptSimulationParameter(sim_handle_self,'forceVectorScaling') forceVectorWidth=simGetScriptSimulationParameter(sim_handle_self,'forceVectorWidth') contactPointSize=simGetScriptSimulationParameter(sim_handle_self,'contactPointSize') overlayDisplay=simGetScriptSimulationParameter(sim_handle_self,'overlayDisplay') if (overlayDisplay) then options=options+sim_drawing_overlay end -- Add a line and a sphere container: lineContainer=simAddDrawingObject(sim_drawing_lines+options,forceVectorWidth,0,-1,1000,black,black,black,purple) sphereContainer=simAddDrawingObject(sim_drawing_spherepoints+options,contactPointSize,0,-1,1000,black,black,black,lightBlue) ObjectHandle = simGetObjectHandle("cube") end if (sim_call_type==sim_childscriptcall_cleanup) then -- Remove the containers: simRemoveDrawingObject(lineContainer) simRemoveDrawingObject(sphereContainer) end if (sim_call_type==sim_childscriptcall_sensing) then -- empty the containers: simAddDrawingObjectItem(lineContainer,nil) simAddDrawingObjectItem(sphereContainer,nil) -- Fill the containers with contact information: index=0 sum = 0 while (true) do objectsInContact,contactPt,forceDirectionAndAmplitude=simGetContactInfo(sim_handle_all,ObjectHandle,index) if (objectsInContact) then sum = sum + math.sqrt(forceDirectionAndAmplitude[1]*forceDirectionAndAmplitude[1]+forceDirectionAndAmplitude[2]*forceDirectionAndAmplitude[2]+forceDirectionAndAmplitude[3]*forceDirectionAndAmplitude[3]) if (index == 39) then str = string.format("%f N", sum/10) simAddStatusbarMessage(str) -- Adds a message to the status bar end line={contactPt[1],contactPt[2],contactPt[3],0,0,0} line[4]=contactPt[1]+forceDirectionAndAmplitude[1]*forceVectorScaling line[5]=contactPt[2]+forceDirectionAndAmplitude[2]*forceVectorScaling line[6]=contactPt[3]+forceDirectionAndAmplitude[3]*forceVectorScaling simAddDrawingObjectItem(lineContainer,line) line[4]=contactPt[1]-forceDirectionAndAmplitude[1]*forceVectorScaling line[5]=contactPt[2]-forceDirectionAndAmplitude[2]*forceVectorScaling line[6]=contactPt[3]-forceDirectionAndAmplitude[3]*forceVectorScaling simAddDrawingObjectItem(lineContainer,line) simAddDrawingObjectItem(sphereContainer,line) index=index+1 else break end end end
再看看接觸力動態變化的例子。在場景中添加一個球體，拖離地面，然后雙擊其圖標打開Dynamic屬性對話框。點擊"Edit material"按鈕，在對應的物理引擎下修改材料屬性。仿真時使用的是Bullet V2.78，這里修改了恢復參數（值越大材料彈性越大）和摩擦參數。注意，當計算兩個物體之間的摩擦時，必須聯合二者的摩擦參數（Two colliding objects will have a combined friction value of value1*value2. This does not correspond to the real friction coefficient）。恢復值通常設置在0到1之間，0意味着小球落地不會彈起，即非彈性碰撞（塑性碰撞，碰撞結束時物體變形無任何恢復，動能全部消耗在碰撞過程中），1為完全彈性碰撞。恢復系數一般由實驗方法確定，工程中材料的恢復系數可以在工程手冊中查到。（Higher restitution values tend to make collisions appear elastic. This does not correspond to the real restitution coefficient. ）
拖入contact display模型顯示接觸力，開始仿真。小球從一定高度自由落下，碰到地面后反彈幾次，最后靜止，下面的動圖可以看出接觸力的變化：
下面一個例子中機械臂末端執行器上裝有力傳感器，用於測量機械手托舉重物的重量。兩個關節采用PID控制使機械臂保持在水平位置，單擊自定義界面上的Add weight按鈕會添加新的重物，對應傳感器測得的重量和關節保持力矩都會增大。另外也直接使用了simGetContactInfo函數獲取接觸力信息並在狀態欄輸出，從下圖可以看到隨着重物的增加接觸力逐漸變大：
Custom UI:
function buttonClick(ui, id) NewObjectHandle = simCreatePureShape(0, 10 , {0.1,0.1,0.1}, 1) simSetObjectPosition(NewObjectHandle, LastHandle, {0,0,0.2}) LastHandle = NewObjectHandle end if (sim_call_type==sim_childscriptcall_initialization) then InitialObjectHandle = simGetObjectHandle('Cuboid') LastHandle = InitialObjectHandle xml = [[ <ui> <button text="Add weight" onclick='buttonClick'/> </ui> ]] ui=simExtCustomUI_create(xml) end if (sim_call_type==sim_childscriptcall_actuation) then end if (sim_call_type==sim_childscriptcall_sensing) then end if (sim_call_type==sim_childscriptcall_cleanup) then simExtCustomUI_destroy(ui) end
ContactInfoDIsplay:
if (sim_call_type==sim_childscriptcall_initialization) then black={0,0,0} purple={1,0,1} lightBlue={0,1,1} options=0 forceVectorScaling=simGetScriptSimulationParameter(sim_handle_self,'forceVectorScaling') forceVectorWidth=simGetScriptSimulationParameter(sim_handle_self,'forceVectorWidth') contactPointSize=simGetScriptSimulationParameter(sim_handle_self,'contactPointSize') overlayDisplay=simGetScriptSimulationParameter(sim_handle_self,'overlayDisplay') if (overlayDisplay) then options=options+sim_drawing_overlay end -- Add a line and a sphere container: lineContainer=simAddDrawingObject(sim_drawing_lines+options,forceVectorWidth,0,-1,1000,black,black,black,purple) sphereContainer=simAddDrawingObject(sim_drawing_spherepoints+options,contactPointSize,0,-1,1000,black,black,black,lightBlue) ObjectHandle = simGetObjectHandle("Pad") end if (sim_call_type==sim_childscriptcall_cleanup) then -- Remove the containers: simRemoveDrawingObject(lineContainer) simRemoveDrawingObject(sphereContainer) end if (sim_call_type==sim_childscriptcall_sensing) then -- empty the containers: simAddDrawingObjectItem(lineContainer,nil) simAddDrawingObjectItem(sphereContainer,nil) -- Fill the containers with contact information: index=0 sum = 0 while (true) do objectsInContact,contactPt,forceDirectionAndAmplitude=simGetContactInfo(sim_handle_all,ObjectHandle,index) if (objectsInContact) then sum = sum + forceDirectionAndAmplitude[3] if (index == 39) then str = string.format("%f N", sum/10) simAddStatusbarMessage(str) -- Adds a message to the status bar end line={contactPt[1],contactPt[2],contactPt[3],0,0,0} line[4]=contactPt[1]-forceDirectionAndAmplitude[1]*forceVectorScaling line[5]=contactPt[2]-forceDirectionAndAmplitude[2]*forceVectorScaling line[6]=contactPt[3]-forceDirectionAndAmplitude[3]*forceVectorScaling simAddDrawingObjectItem(lineContainer,line) simAddDrawingObjectItem(sphereContainer,line) index=index+1 else break end end end
參考：
Video Game Physics Tutorial - Part III: Constrained Rigid Body Simulation