歐惠棠 李晉芳 莫建清

(廣東工業(yè)大學(xué)機(jī)電工程學(xué)院 廣州 510006)

引言

虛擬醫(yī)療是應(yīng)用虛擬現(xiàn)實(shí)(virutalreatility,VR)技術(shù),構(gòu)建虛擬的人體模型器官以及手術(shù)，提高虛擬實(shí)境的真實(shí)感。借助于虛擬外設(shè)，可以使人們更逼真地學(xué)習(xí)醫(yī)療知識(shí)以及治病救人的虛擬現(xiàn)實(shí)應(yīng)用，在未來(lái)有著非常樂(lè)觀的前景。其中，手術(shù)訓(xùn)練與虛擬現(xiàn)實(shí)相結(jié)合的虛擬手術(shù)是現(xiàn)代手術(shù)訓(xùn)練的新趨勢(shì)。虛擬手術(shù)仿真利用醫(yī)學(xué)生物信息和計(jì)算機(jī)仿真技術(shù)以及其他仿真模擬器械，為醫(yī)療工作者提供一個(gè)虛擬且逼真的手術(shù)場(chǎng)景，使其借助這些虛擬場(chǎng)景進(jìn)行手術(shù)計(jì)劃訓(xùn)練和研究，積累經(jīng)驗(yàn)，掌握技巧，以便在實(shí)際的手術(shù)過(guò)程中完成復(fù)雜的手術(shù)操作，為病患減少手術(shù)后的痛苦和手術(shù)風(fēng)險(xiǎn)[1]。通過(guò)這一手段，無(wú)疑可使訓(xùn)練更真實(shí)、準(zhǔn)確、可靠。

國(guó)內(nèi)外對(duì)于虛擬手術(shù)仿真中碰撞力的研究，著重于變形體形變技術(shù)的實(shí)現(xiàn)上，如Montgomery等提出的現(xiàn)在非常流行的彈簧-質(zhì)點(diǎn)網(wǎng)絡(luò)(spring-mass networks)[2]，F(xiàn)elippa等提出的基于共旋轉(zhuǎn)模型的有限元方法(FEM)[3]，吳涓等的一種基于沿徑向方向分割為呈同心圓分布的彈簧-質(zhì)點(diǎn)模型及實(shí)時(shí)力學(xué)響應(yīng)算法[4]，陳衛(wèi)東等提出的虛擬體彈簧的變形算法[5]等。這些研究?jī)H僅著重于變形體內(nèi)部力的計(jì)算上，并沒(méi)有系統(tǒng)地給出碰撞檢測(cè)、碰撞力的計(jì)算、時(shí)間積分等整個(gè)仿真流程。國(guó)內(nèi)的研究?jī)H僅給出當(dāng)模型發(fā)生形變時(shí)內(nèi)力與粒子群位移之間的關(guān)系，卻沒(méi)有給出與手術(shù)器械模型發(fā)生碰撞時(shí)如何計(jì)算出碰撞力并使變形體發(fā)生變形。國(guó)外Hasegawa等提出懲罰方法來(lái)求解碰撞接觸力[6]，但會(huì)導(dǎo)致一系列不穩(wěn)定性問(wèn)題。Anitescu等提出k邊形金字塔算法[7]，但時(shí)間復(fù)雜度較大，導(dǎo)致運(yùn)算時(shí)間過(guò)長(zhǎng)。

本研究提出虛擬手術(shù)訓(xùn)練中帶摩擦的碰撞力計(jì)算及其仿真方法，闡述粒子系統(tǒng)的建立、碰撞檢測(cè)得到碰撞力方向、隱式歐拉積分、碰撞約束方程(線性互補(bǔ)問(wèn)題)、摩擦約束方程的構(gòu)建、求解方程得出碰撞力大小這一系列的完整流程。該方法適用于剛體和變形體(彈簧質(zhì)點(diǎn)網(wǎng)絡(luò)和有限元)物理模型仿真，不僅實(shí)現(xiàn)了碰撞功能，而且還添加了摩擦特性，能夠在一定接觸點(diǎn)數(shù)目的場(chǎng)景下滿足交互的真實(shí)性、穩(wěn)定性和實(shí)時(shí)性的要求。

在計(jì)算機(jī)圖形學(xué)中，物理仿真通過(guò)粒子系統(tǒng)[8]實(shí)現(xiàn)。粒子系統(tǒng)可以建立剛體和變形體等物理對(duì)象，并能實(shí)現(xiàn)其物理行為。本文所述的碰撞研究對(duì)剛體和變形體均適用。對(duì)于碰撞物理仿真的流程，已有相關(guān)的研究。下面通過(guò)約束動(dòng)力學(xué)構(gòu)建碰撞與摩擦相關(guān)約束方程并求解，然后通過(guò)實(shí)驗(yàn)驗(yàn)證是否計(jì)算正確和性能是否符合要求。

1 材料和方法

1.1 材料

采用自行設(shè)計(jì)的虛擬手術(shù)仿真系統(tǒng)，加入帶摩擦的碰撞力計(jì)算算法，進(jìn)行藍(lán)鉗模型和半月板模型的碰撞測(cè)試。實(shí)驗(yàn)環(huán)境為Win10 64 位系統(tǒng)，處理器為Intel(R) Core(TM) i5-3470CPU @ 3.20 GHz，顯卡為NVIDIA Quadro 600。先分別通過(guò)藍(lán)鉗與半月板模型的單接觸面和雙接觸面碰撞進(jìn)行10次實(shí)驗(yàn)，驗(yàn)證所計(jì)算出的碰撞力方向和響應(yīng)時(shí)間是否達(dá)到要求，再驗(yàn)證100次測(cè)試實(shí)驗(yàn)中，單接觸面接觸碰撞從1個(gè)接觸點(diǎn)到45個(gè)接觸點(diǎn)同時(shí)碰撞的平均碰撞時(shí)間。

1.2 方法

計(jì)算模型的碰撞接觸力的基本仿真流程如圖1所示。

圖1 碰撞仿真流程Fig.1 Schematic diagram of collision simulation

為了計(jì)算接觸力，必須先通過(guò)碰撞檢測(cè)，檢測(cè)出潛在的碰撞對(duì)象，當(dāng)模型發(fā)生碰撞后，精確計(jì)算出碰撞點(diǎn)的位置和接觸力的方向(大小未知)，然后構(gòu)建相應(yīng)的約束方程，再利用相應(yīng)的數(shù)值方法進(jìn)行求解。求解后，通過(guò)圖形界面驗(yàn)證計(jì)算出的反饋力的方向和大小是否符合要求，通過(guò)20次的碰撞測(cè)試，并作出碰撞力的曲線，驗(yàn)證接觸力的響應(yīng)時(shí)間和穩(wěn)定性是否符合要求。通過(guò)細(xì)分半月板模型的網(wǎng)格并通過(guò)100次的實(shí)驗(yàn)，驗(yàn)證單接觸面接觸從1個(gè)接觸點(diǎn)到50個(gè)接觸點(diǎn)同時(shí)接觸的平均碰撞時(shí)間。

1.2.1碰撞檢測(cè)

如圖2所示，通過(guò)碰撞檢測(cè)，可以得出兩個(gè)即將碰撞的模型的潛在碰撞點(diǎn)P和Q?？梢匀我膺x取法向量，本研究選取模型D1上P點(diǎn)對(duì)應(yīng)的朝內(nèi)的法向量，記為n。

圖2 兩個(gè)物體之間的碰撞Fig.2 The collision between two objects

定義兩個(gè)模型在P點(diǎn)上的空隙為

(1)

fn21(P)+fn12(Q)=0

(2)

由Signorini法則[9]，對(duì)于碰撞行為，δn(P)與fn21(P)非零并且具有互補(bǔ)關(guān)系，記為

0≤δn(P)⊥fn21(P)≥0

(3)

互補(bǔ)關(guān)系是指δn(P)與fn21(P)中必有一個(gè)為零[10]。

1.2.2插值

在計(jì)算機(jī)圖形系統(tǒng)中，模型表面為三角面片構(gòu)成的網(wǎng)格。模型之間的碰撞化歸為兩個(gè)三角形之間的標(biāo)準(zhǔn)碰撞，如圖3所示。碰撞檢測(cè)出來(lái)的碰撞點(diǎn)可能位于三角形的內(nèi)部或邊線上，但由于在粒子系統(tǒng)中，粒子通常是三角面片的頂點(diǎn)，因此必須把碰撞點(diǎn)的碰撞力等效到三角形的頂點(diǎn)上。

圖3 兩個(gè)三角形碰撞的全部5種情況Fig.3 All five cases of the collision between two triangles

根據(jù)解析幾何相關(guān)知識(shí)，三角形所在平面的任意點(diǎn)都能表示為頂點(diǎn)的加權(quán)平均值，其中的權(quán)稱為重心坐標(biāo)。在圖4中，對(duì)于點(diǎn)P，有

P=αA1+βB1+γC1

(4)

式中，α、β、γ為重心坐標(biāo)分量。

式(4)中兩邊取微位移，有

u(P)=αu(A1)+βu(B1)+γu(C1)

(5)

圖4 碰撞點(diǎn)在三角形內(nèi)部Fig.4 The contact point located in the interior of atriangle

由虛功原理，有

uT(P)fP=uT(A1)fA1+uT(B1)fB1+uT(C1)fC1

(6)

從而

uT(A1)(fA1-αfP)+uT(B1)(fB1-βfP)+

uT(C1)(fC1-γfP)=0

(7)

得出

fA1=αfPfB1=βfPfC1=γfP

(8)

1.2.3構(gòu)建線性互補(bǔ)問(wèn)題

在圖2中，設(shè)D1、D2中所有粒子的位置分別為q1、q2，三角形碰撞面片的頂點(diǎn)與q1、q2有一定的映射關(guān)系，稱為接觸蒙皮(contact skinning)，不詳細(xì)討論。接觸蒙皮能使得兩個(gè)模型之間三角面片頂點(diǎn)的約束轉(zhuǎn)化為模型之間粒子群的約束。因此，δn(P)≥0可表示為

Ψ(q1,q2)≥0

(9)

稱為單面約束，Ψ僅與位置有關(guān)，稱為完整約束[11]。

粒子系統(tǒng)的物理仿真實(shí)質(zhì)上是求解牛頓第二定律常微分方程(ODE)有

(10)

式中，fext為外部力，fint為內(nèi)部力，fcon為約束力。

令fext+fint=f(q,v)，fcon=HTλ，其中q為粒子群的位置，v為粒子群的速度，H為約束力的方向，λ為約束力的大小。

使用隱式歐拉方法來(lái)求解常微分方程，有

vt+h-vt=hat+h

(11)

qt+h-qt=hvt+h

(12)

(13)

式中，t為上一個(gè)時(shí)間步的時(shí)刻，h為隱式歐拉方法時(shí)間步的步長(zhǎng)。

(M-hB-h2K)dv=hf(xt,vt)+h2Kvt+HTλh

(14)

式中，dv=vt+h-vt。

令A(yù)=M-hB-h2K,b=hf(qt,vt)+h2Kvt，新的λ=λh。則有Adv=b+HTλ，從而對(duì)兩個(gè)粒子群q1、q2之間的約束，有

(15)

若取λ=0，求得的解為沒(méi)有約束的運(yùn)動(dòng)狀態(tài)，稱為自由運(yùn)動(dòng)解，為dvfree1和dvfree2。

事實(shí)上，dv=dvfree+dvcor，因此，

(16)

最后得到

(17)

對(duì)于線性互補(bǔ)問(wèn)題，其求解過(guò)程可由圖5直觀表示。

1.2.4建立摩擦約束

圖5 線性互補(bǔ)問(wèn)題求解。(a)未發(fā)生碰撞；(b)發(fā)生碰撞Fig.5 Linear complementar yproblem solving.(a)Collision does not occur； (b)Collision occurs

對(duì)于摩擦約束，僅知道摩擦力位于所在碰撞點(diǎn)的切平面上，而在切平面的哪個(gè)方向未知，需要求解。令摩擦力為ft，滑動(dòng)摩擦位移為δt，是一個(gè)二維向量。

碰撞約束方程與摩擦約束方程合并，得到總約束方程為

(18)

其中，接觸約束方程為線性互補(bǔ)問(wèn)題，有

0≤fn⊥δn≥0

(19)

對(duì)于摩擦約束方程，要滿足干摩擦法則(Coulomb’s law)，有

對(duì)于摩擦約束，其求解可由圖6直觀表示。

圖6 摩擦約束求解。(a)靜摩擦；(b)動(dòng)摩擦Fig.6 Friction constrain tsolving.(a) static friction； (b) dynamical friction

1.2.5建立總約束方程

對(duì)于接觸面約束，總的約束方程為式(18)。對(duì)于復(fù)雜的手術(shù)交互，除了單接觸面碰撞之外，還會(huì)涉及雙接觸面碰撞，如軟組織夾取交互。若雙接觸面碰撞仍按照單接觸面碰撞的情況處理，則一個(gè)接觸面碰撞響應(yīng)會(huì)影響另一個(gè)接觸面的碰撞響應(yīng)，導(dǎo)致兩個(gè)接觸面不能穩(wěn)定地產(chǎn)生碰撞壓力。靈活地改變約束方程，能使約束方程仍然適用于雙接觸面的情形。

如果碰撞檢測(cè)出軟組織的同時(shí)與兩個(gè)接觸面發(fā)生碰撞，則約束方程發(fā)生改變。如圖7所示，圖中三角面片代表接觸區(qū)域。(a)中軟組織模型與夾子模型的一個(gè)面接觸，此時(shí)采用單接觸面約束求解方法。(b)中軟組織模型與夾子模型的兩個(gè)面同時(shí)接觸，則原來(lái)的單接觸面狀態(tài)轉(zhuǎn)變?yōu)殡p接觸面狀態(tài)，使用雙接觸面約束求解方法。

圖7 接觸狀態(tài)的轉(zhuǎn)變。(a)單接觸面接觸；(b)雙接觸面接觸Fig.7 Change of contact status.(a)Single contact surface collision； (b)Double contact surface collision

軟組織模型碰撞點(diǎn)所在的三角面片頂點(diǎn)與夾子碰撞面建立新的線性互補(bǔ)，有

0≤δn⊥(ξfn+∑ξadjfnadj)≥0

(22)

建立基于頂點(diǎn)領(lǐng)域的摩擦方程，有

δt=0?‖ft‖<μ‖ξfn+∑ξadjfnadj‖

(23)

(24)

式中，fnadj為頂點(diǎn)的鄰接頂點(diǎn)的法向約束力，ξ為權(quán)重。

當(dāng)fn≤0，ξ=0；當(dāng)fn>0，ξ=fn/N(N為頂點(diǎn)及其鄰接頂點(diǎn)的數(shù)目)。

由于新的約束方程允許fn≤0，并且把三角面片頂點(diǎn)及其鄰接頂點(diǎn)法向力的綜合作用考慮進(jìn)來(lái)，可以有效提高交互的穩(wěn)定性，并且實(shí)質(zhì)上更符合于現(xiàn)實(shí)的接觸情況。

從而雙接觸面接觸總的約束方程為

(25)

式中，m為接觸的三角面片頂點(diǎn)的數(shù)目。

1.2.6求解總約束方程

對(duì)于式(25)，使用類高斯-賽德?tīng)査惴?Gauss-Seidel-like algorithm)求解。在粒子仿真計(jì)算力學(xué)領(lǐng)域[13]中，該算法能保證較好的收斂性。先對(duì)接觸逐個(gè)求解，然后不停地迭代，直到數(shù)值收斂。對(duì)于m個(gè)瞬時(shí)摩擦接觸，先考慮其中一個(gè)摩擦接觸α，對(duì)其進(jìn)行求解，有

(26)

求出fα，其他接觸也以此進(jìn)行計(jì)算，然后不斷迭代求出新的fα，直到數(shù)值收斂。

(27)

本研究?jī)H給出雙接觸面接觸約束求解算法，單接觸面接觸求解算法與雙接觸面接觸類似，并且要簡(jiǎn)單，由于篇幅所限，這是不再贅述。由于雙接觸面接觸要求查詢?nèi)敲嫫旤c(diǎn)的鄰接頂點(diǎn)，故模型的網(wǎng)格要求為流形網(wǎng)格，并采用半邊數(shù)據(jù)結(jié)構(gòu)，可以使查詢操作的時(shí)間復(fù)雜度為O(1)。算法如下：

Input:(δfree)3 m×1,[W]3 m×3m

Output:(f)3 m×1

設(shè)置閾值ε1、ε2、ε3

k=0

fori=1…mdo

(fi(0))(3×1)=0

Λmin/max=eig([Wii]tt) 每個(gè)([Wii]tt)均有兩個(gè)特征值

end

repeat

k=k+1

foreachi=1…mdo

(δ)(3×1)=(δfree i)(3×1)

end

foreachj=1…i-1do

(δ)(3×1)+=[Wij]3×3(fj(k))3×1

end

foreachj=i…mdo

(δ)(3×1)+=[Wij](3×3)(fj(k-1))(3×1)

end

if(fi)nadj>ε1then

(fi(k))t=(fi(k-1))t-(δ)t/Λi

else

(fi(k))=0

end

foreachi=1…mdo

if(fi)n>ε3thenξ=1

elseξ=0

end

(fi)nadj=ξ(fi)n

((fik)nadj)l×1=adjacent Vertices Normal Forces(i)

foreachs=1…ldo

if(fis)nadj>ε3thenξs=1

elseξs=0

end

(fi)nadj+=ξs(fis)nadj

end

(fi)nadj/=l+1

if(fi)nadj≤ε3

(fi(k))t=0

continue

end

if‖(fi(k))t‖>μ(fi)nadjthen

(fi(k))t×=μ(fi)nadj/‖(fi(k))t‖

end

end

1.2.7實(shí)驗(yàn)方法

在碰撞試驗(yàn)中，類高斯-賽德?tīng)査惴ㄊ諗块撝郸?、ε2、ε3取0.001，摩擦系數(shù)μ取0.5。半月板模型設(shè)定楊氏模量為59 MPa，泊松比為0.45[14]。實(shí)驗(yàn)方案為先分別進(jìn)行10次單接觸面和雙接觸面接觸碰撞實(shí)驗(yàn)。單接觸面接觸時(shí)，藍(lán)鉗模型簡(jiǎn)單地碰撞半月板模型；雙接觸面接觸時(shí)，藍(lán)鉗模型夾取碰撞模型。碰撞時(shí)，分析碰撞力的大小和方向、響應(yīng)時(shí)間、穩(wěn)定時(shí)候的起伏是否符合要求。然后，再細(xì)分半月板模型的網(wǎng)格，使得碰撞時(shí)的接觸點(diǎn)增加，進(jìn)行100次實(shí)驗(yàn)，計(jì)算1個(gè)接觸點(diǎn)到45個(gè)接觸點(diǎn)同時(shí)接觸的平均運(yùn)算時(shí)間，其運(yùn)算時(shí)間的記錄方法如圖8所示。

圖8 碰撞時(shí)間測(cè)試流程Fig.8 Diagram of the collision time test

1.2.8統(tǒng)計(jì)學(xué)分析

單接觸面接觸中，以膝關(guān)節(jié)模型的正方向(髕骨方向)為基準(zhǔn)，分別在豎直和水平方向從-90°、45°、0°、45°、90°各進(jìn)行1次碰撞實(shí)驗(yàn)，計(jì)算平均的響應(yīng)時(shí)間和穩(wěn)定時(shí)的起伏。在雙接觸面接觸中，從水平方向-90°、45°、0°、45°、90°進(jìn)行10次夾取試驗(yàn)，計(jì)算平均的響應(yīng)時(shí)間和穩(wěn)定時(shí)的起伏。在多接觸點(diǎn)平均碰撞力運(yùn)算時(shí)間的測(cè)試中，為了能使足夠的接觸點(diǎn)接觸，接觸的方向任意選取，然后計(jì)算平均計(jì)算時(shí)間。

2 結(jié)果

圖9為筆者構(gòu)建的膝關(guān)節(jié)鏡手術(shù)場(chǎng)景，在場(chǎng)景中，藍(lán)鉗模型與半月板模型發(fā)生碰撞。碰撞接觸力通過(guò)上述方法得到計(jì)算并繪制，由圖可知接觸力的方向符合實(shí)際情況。

圖9 場(chǎng)景中碰撞力的繪制Fig.9 Contact force rendering in the scene

圖10為10次單接觸面碰撞實(shí)驗(yàn)中其中一次的碰撞力的曲線。由圖可知，模型在1和2 s的時(shí)候發(fā)生了碰撞，x、y、z方向力的大小變化趨勢(shì)無(wú)顯著差異，開(kāi)始碰撞時(shí)，接觸力急劇上升，隨著半月板模型的變形行為趨于穩(wěn)定，碰撞力也趨于穩(wěn)定。在2 s的時(shí)候，藍(lán)鉗模型在z軸負(fù)方向前進(jìn)0.2個(gè)單位長(zhǎng)度，穩(wěn)定后碰撞力增加。在10次的實(shí)驗(yàn)中，模型碰撞的平均響應(yīng)時(shí)間為0.02 s，接觸力穩(wěn)定后的平均起伏范圍為±0.02個(gè)單位。

圖10 碰撞力的分析Fig.10 Analysis of the contact force

圖11 通過(guò)摩擦力夾取變形體Fig.11 Grasping a deformable object by friction

圖11為半月板模型夾取的效果。在不同的時(shí)間段，其狀態(tài)從左到右、從上到下分別為：準(zhǔn)備夾取階段、夾取階段、向右拖動(dòng)階段、向下拖動(dòng)階段、釋放階段。

圖12為10次夾取實(shí)驗(yàn)中的一次藍(lán)鉗下顎碰撞力曲線。在6.12 s的時(shí)候，藍(lán)鉗夾取半月板，接觸力在y方向發(fā)生突變，隨后趨于穩(wěn)定；在8.78 s的時(shí)候，藍(lán)鉗拉動(dòng)半月板運(yùn)動(dòng)，摩擦使接觸力在z軸方向上發(fā)生突變；在10.78 s的時(shí)候，藍(lán)鉗釋放半月板，接觸力回復(fù)為零。在10次的實(shí)驗(yàn)中，模型碰撞的平均響應(yīng)時(shí)間為0.02 s，接觸力穩(wěn)定后的平均起伏范圍為±0.02個(gè)單位。

圖13為100次單接觸面碰撞下1～45個(gè)接觸點(diǎn)同時(shí)接觸時(shí)的平均運(yùn)算時(shí)間散點(diǎn)。運(yùn)算時(shí)間的測(cè)試結(jié)果顯示，1個(gè)接觸的運(yùn)算時(shí)間約為0.9 ms，11個(gè)接觸點(diǎn)時(shí)平均每個(gè)接觸點(diǎn)的運(yùn)算時(shí)間為1.9 ms。隨著接觸點(diǎn)數(shù)目的增長(zhǎng)，運(yùn)算時(shí)間也急劇上升。到44個(gè)接觸點(diǎn)時(shí)，平均每個(gè)接觸的運(yùn)算時(shí)間為8 ms，實(shí)時(shí)性有所降低。

圖13 不同接觸點(diǎn)數(shù)目的運(yùn)算時(shí)間Fig.13 Computation time of various numbers of contact points

3 討論

單接觸面接觸的實(shí)驗(yàn)結(jié)果表明，單接觸面接觸過(guò)程的運(yùn)算速度快，計(jì)算出來(lái)的接觸力穩(wěn)定，能很好地達(dá)到虛擬手術(shù)中的實(shí)時(shí)性和真實(shí)性要求。

雙接觸面接觸的實(shí)驗(yàn)結(jié)果表明，在雙接觸面的情況下能計(jì)算出接觸力和摩擦力，穩(wěn)定性能得到保證，得到良好的摩擦效果，仍能滿足虛擬手術(shù)中的實(shí)時(shí)性和真實(shí)性要求。

在多接觸點(diǎn)同時(shí)接觸的情形下，碰撞力的計(jì)算時(shí)間與接觸點(diǎn)的數(shù)量有直接的關(guān)系，會(huì)隨著接觸點(diǎn)的增加而增長(zhǎng)。在多模型的場(chǎng)景中，碰撞的發(fā)生會(huì)較多，從而碰撞點(diǎn)增多，碰撞的時(shí)間也越長(zhǎng)。但是，若模型不發(fā)生碰撞，那么模型數(shù)量影響的僅僅是下一個(gè)時(shí)間步粒子群位置的計(jì)算時(shí)間，接觸計(jì)算不會(huì)發(fā)生。

在本次多接觸點(diǎn)接觸實(shí)驗(yàn)中，采用類高斯-賽德?tīng)査惴ㄇ蠼饩€性互補(bǔ)問(wèn)題。類高斯-賽德?tīng)査惴ǖ臅r(shí)間復(fù)雜度為O(m2)，m為接觸點(diǎn)的數(shù)目，顯然求解的時(shí)間隨著接觸點(diǎn)的增加而增加，與本次實(shí)驗(yàn)的結(jié)果基本吻合。這表明，隨著接觸點(diǎn)的增加，仍然會(huì)使總的求解時(shí)間上升，從而影響實(shí)時(shí)性。國(guó)外有一種求解帶摩擦接觸的k邊形金字塔算法[7]，其時(shí)間復(fù)雜度為O(k×m2)，m為接觸點(diǎn)的數(shù)目，k為近似摩擦圓錐底面的多邊形邊數(shù)。顯然，類高斯-賽德?tīng)査惴ǜ邇?yōu)勢(shì)。

本研究表明，手術(shù)交互的方式復(fù)雜而多樣，即便是摩擦接觸交互，對(duì)不同的接觸方式也不能使用單一不變的方法。本研究使用了基于約束的方法來(lái)處理摩擦接觸交互問(wèn)題：對(duì)于單接觸面接觸，使用單面約束和摩擦約束，簡(jiǎn)單而有效；對(duì)于雙接觸面接觸，則原本的約束方程不再適用，改為基于三角面片頂點(diǎn)及其鄰接頂點(diǎn)的單面約束和摩擦約束。雙接觸面接觸的情形比單接觸面的情況要復(fù)雜一些，但能提高交互的穩(wěn)定性，也沒(méi)有跳出約束方程這一方法的框架。由此表明，使用基于約束方程的方法，能夠應(yīng)對(duì)復(fù)雜多變的交互，即使交互的形式不同，仍能通過(guò)適當(dāng)?shù)馗淖兗s束方程來(lái)實(shí)現(xiàn)。

4 結(jié)論

本研究提出了虛擬手術(shù)中帶摩擦的接觸力的計(jì)算方法。首先，建立計(jì)算機(jī)圖形學(xué)物理仿真中最為常用的粒子系統(tǒng)，通過(guò)碰撞檢測(cè)出潛在的碰撞對(duì)象，精確計(jì)算碰撞接觸點(diǎn)和接觸力的方向，總體的常微分方程通過(guò)隱式歐拉方法進(jìn)行時(shí)間積分；其次，建立碰撞接觸和摩擦接觸的約束方程，并通過(guò)類高斯-賽德?tīng)柗椒ㄇ蠼?；再次，分別進(jìn)行10次單接觸面和雙接觸面接觸碰撞實(shí)驗(yàn)，碰撞時(shí)分析碰撞力的大小和方向、響應(yīng)時(shí)間、穩(wěn)定時(shí)候的起伏是否符合要求；最后，再細(xì)分半月板模型的網(wǎng)格，使得碰撞時(shí)的接觸點(diǎn)增加，進(jìn)行100次實(shí)驗(yàn)，計(jì)算1～45個(gè)接觸點(diǎn)同時(shí)接觸的平均運(yùn)算時(shí)間。

實(shí)驗(yàn)結(jié)果表明，通過(guò)這一系列方法，在建立的變形體夾取場(chǎng)景中，能繪制出符合實(shí)際情況的力的大小和方向。該方法在單接觸面和雙接觸面接觸的情況下，滿足虛擬手術(shù)的真實(shí)性和實(shí)時(shí)性要求。在多接觸點(diǎn)接觸的情況下，隨著接觸點(diǎn)數(shù)目的增長(zhǎng)，運(yùn)算時(shí)間也急劇上升，然而在11個(gè)接觸點(diǎn)的情況下，仍能滿足實(shí)時(shí)性要求。求解約束方程所用的高斯-賽德?tīng)柗椒ň哂懈偷臅r(shí)間復(fù)雜度O(m2)，而且有很好的靈活性，通過(guò)稍加修改即可用于雙接觸面接觸的情形，并且能很好地保持其性能。不過(guò)，在接觸點(diǎn)的數(shù)目上升到一定程度后，還是會(huì)對(duì)仿真的實(shí)時(shí)性產(chǎn)生影響。

未來(lái)將會(huì)與自行設(shè)計(jì)的力反饋設(shè)備進(jìn)行虛耦合，建立一整套力反饋交互系統(tǒng)。