本發(fā)明涉及一種人體組織有限元仿真技術(shù)，特別涉及一種軟組織有限元模型的加速計算方法領(lǐng)域。

背景技術(shù)：

虛擬手術(shù)培訓(xùn)系統(tǒng)能夠提高醫(yī)生的手術(shù)技能，從而降低手術(shù)失敗的風(fēng)險。構(gòu)建一個可靠的虛擬手術(shù)培訓(xùn)系統(tǒng)，不僅需要精確地反映人體組織真實的力學(xué)特性，還有保證其仿真的實時性。人體軟組織具有復(fù)雜的力學(xué)特性，而利用有限元法能夠精確地預(yù)測人體組織的變形，但由于其較長的計算時間，難以在虛擬手術(shù)培訓(xùn)系統(tǒng)中直接應(yīng)用。

對現(xiàn)有技術(shù)的文獻(xiàn)檢索發(fā)現(xiàn)，閆楨南在上海交通大學(xué)的碩士學(xué)位論文《軟組織非線性形變快速模型與應(yīng)用》中，基于有限元法將軟組織模型分為操作區(qū)域和非操作區(qū)域，并只對操作區(qū)域的形變進(jìn)行計算的方法來減少軟組織模型計算規(guī)模，從而達(dá)到加速仿真的效果。而這種方法需要用戶在模擬手術(shù)演練前，根據(jù)自己以往的經(jīng)驗，對軟組織的模型進(jìn)行手動劃分操作區(qū)域，因此缺乏靈活性，尤其是對手術(shù)路徑較長的情況，其加速效果并不理想。

綜上所述，目前在建立虛擬手術(shù)培訓(xùn)系統(tǒng)的過程中，軟組織的實時仿真仍然是一項巨大的挑戰(zhàn)。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的在于為了克服現(xiàn)有技術(shù)中的不足，提供一種軟組織有限元模型的加速計算方法，能夠同時實現(xiàn)虛擬手術(shù)培訓(xùn)系統(tǒng)的真實性和實時性。

本發(fā)明為解決上述技術(shù)問題采取的技術(shù)方案是：

一種軟組織有限元模型的加速計算方法，其特征在于包括以下步驟：

步驟1）：根據(jù)軟組織的網(wǎng)格模型，創(chuàng)建數(shù)組A，將軟組織模型的每個四面體單元序號與每個四面體單元所包含的四個三角形單元序號相連，使得任意一個四面體單元序號可得出其包含的四個三角形單元序號。同時創(chuàng)建數(shù)組B，將軟組織模型的每個三角形單元序號與共用每個三角形單元的兩個相鄰四面體單元的序號相連，使得通過任意一個三角形單元序號可得出共用其三角形單元的兩個相鄰四面體單元序號，如果其中某一三角形單元為軟組織模型的邊界單元，則只與包含當(dāng)前三角形單元的四面體單元序號相連。

步驟2）：篩選參與計算的單元的集合，以受力點所在單元的內(nèi)心為基準(zhǔn)點，給出如下所示的篩選條件的關(guān)系式：

(1)

式中，l為當(dāng)前四面體單元的內(nèi)心與受力點所在單元的內(nèi)心之間的距離，a為篩選系數(shù)，R為受力單元的外接圓半徑；

判斷當(dāng)前的四面體單元是否滿足篩選條件關(guān)系式（1），如果是，將當(dāng)前四面體單元列入?yún)⑴c計算單元的集合，判斷首先從受力點所在四面體單元的相鄰四面體單元開始進(jìn)行，然后再對每個相鄰四面體單元相鄰的四面體單元進(jìn)行判斷，參與計算單元集合的判斷以受力點所在四面體單元為中心，由里到外進(jìn)行，直到不滿足篩選條件關(guān)系式（1）為止。篩選參與計算的單元時，根據(jù)篩選判斷的先后順序，將參與計算的單元集合的四面體單元序號和節(jié)點序號進(jìn)行重新排序。當(dāng)受力點所在四面體單元個數(shù)大于1時，依次進(jìn)行受力點所在四面體單元個數(shù)次的循環(huán)進(jìn)行判斷，并統(tǒng)一封裝到參與計算單元的集合；每個四面體單元的相鄰四面體單元，通過步驟1）中生成的數(shù)組A和數(shù)組B進(jìn)行查找。

步驟3）：創(chuàng)建數(shù)組C，將步驟2）中生成的參與計算單元集合的新四面體單元序號和原模型中的四面體單元序號相連。創(chuàng)建數(shù)組D，將步驟2）中生成的參與計算單元集合的新節(jié)點序號和原模型中的節(jié)點序號相連。

步驟4）：基于步驟2）所生成的參與計算的單元集合，并結(jié)合步驟3）中所生成的數(shù)組C和數(shù)組D，封裝整體剛度矩陣，在施加的外加載荷和邊界條件下，只對軟組織模型中所篩選出的參與計算單元的集合進(jìn)行形變計算。

上述的一種軟組織有限元模型的加速計算方法，所述步驟1）中的數(shù)組A為m行4列的二維數(shù)組，數(shù)據(jù)類型為int型，m為軟組織模型所包含的四面網(wǎng)格單元的個數(shù)，每個行下標(biāo)代表軟組織模型的四面體網(wǎng)格單元的序號，每行的四個列下標(biāo)所對應(yīng)的元素分別存儲每個行下標(biāo)所對應(yīng)四面體單元所包含的四個三角形單元的序號。數(shù)組B為n行2列的二維數(shù)組，數(shù)據(jù)類型為int型，n為軟組織模型所包含的三角形單元的個數(shù)，每個行下標(biāo)代表軟組織模型的三角形網(wǎng)格單元的序號，每行的兩個列下標(biāo)所對應(yīng)的元素分別存儲共用當(dāng)前三角形單元的兩個四面體單元的序號，如果當(dāng)前三角形單元為軟組織模型邊界面上的單元，則第一個元素存儲包含當(dāng)前三角形單元的四面體單元的序號，并且第二個元素為-1。

上述的一種軟組織有限元模型的加速計算方法，所述步驟3)中的數(shù)組C為int型的一維數(shù)組，數(shù)組C的下標(biāo)表示原模型的四面體單元序號，每個參與計算單元的新四面體單元序號分別存儲到對應(yīng)的原模型的四面體單元序號下標(biāo)的數(shù)組C的元素中，而原模型的四面體單元序號中，非參與計算單元下標(biāo)的元素值賦為-1。數(shù)組D為int型的一維數(shù)組，數(shù)組D的下標(biāo)表示原模型的節(jié)點序號，每個參與計算單元的新節(jié)點序號分別存儲到對應(yīng)的原模型節(jié)點序號下標(biāo)的數(shù)組D的元素中，而原模型節(jié)點序號中，非參與計算節(jié)點下標(biāo)的對應(yīng)元素值賦為-1。

上述的一種軟組織有限元模型的加速計算方法，所述步驟4）中軟組織形變計算的過程中，軟組織的模型采用指數(shù)多項式混合形式的超彈性模型，模型的表達(dá)式如下：

(2)

式中，W為應(yīng)變能密度函數(shù)，C₁、β為模型的材料常數(shù)，I₁和I₂為應(yīng)變張量的主不變量；

基于軟組織的指數(shù)多項式混合形式的超彈性模型建立平衡方程，然后施加邊界條件和外加載荷，可進(jìn)行形變計算。在進(jìn)行形變計算的過程中，平衡方程采用完全拉格朗日法，進(jìn)行分步載荷求解，每一載荷步的非線性方程組采用修正的牛頓迭代法進(jìn)行求解。在軟組織有限元計算的過程中，整體剛度矩陣采用GPU計算，非線性方程組采用CPU進(jìn)行求解。

本發(fā)明的有益效果是：

本發(fā)明中，基于指數(shù)多項式混合形式的超彈性模型，進(jìn)行人體軟組織的有限元仿真，能夠在一定程度上保證其仿真精度。然后，根據(jù)圣維南原理，通過給出一個篩選系數(shù)a來篩選軟組織模型中參與計算單元的集合，并只對局部篩選出的參與計算的單元集合進(jìn)行有限元計算，

提高軟組織有限元模型的計算速度。用戶通過合理地選取篩選系數(shù)a的值，能夠有效地減少系統(tǒng)的計算規(guī)模，同時能夠真實地仿真出人體軟組織的生物力學(xué)特性。整個過程中，整體剛度矩陣采用GPU進(jìn)行計算，能夠進(jìn)一步地提高了整個有限元計算的速度。

附圖說明

圖1是本發(fā)明流程框圖，即軟組織加速計算流程圖；

圖2是篩選參與計算單元集合的判斷示意圖；

圖3是本發(fā)明應(yīng)用的平面示意圖。

具體實施方式

具體實施方式一：如圖1所示，本實施方式所述的一種軟組織有限元模型的加速計算方法，包括以下步驟：

步驟1）：根據(jù)軟組織的網(wǎng)格模型中的每個四面體單元序號與三角形單元序號，創(chuàng)建數(shù)組A，將軟組織模型的每個四面體單元序號與每個四面體單元所包含的四個三角形單元序號相連，使得通過軟組織模型的任意一個四面體單元序號可得出其包含的四個三角形單元序號。同時創(chuàng)建數(shù)組B，將軟組織模型的每個三角形單元序號與共用每個三角形單元的兩個相鄰四面體單元序號相連，使得通過軟組織模型的任意一個三角形單元序號可得出共用該三角形單元的兩個相鄰四面體單元序號，如果其中某一三角形單元為軟組織模型的邊界單元，則只與包含當(dāng)前三角形單元的四面體單元序號相連。本步驟在離線狀態(tài)下實現(xiàn)，因此對軟組織計算速度無影響。

步驟2）：篩選參與計算的單元的集合，以受力點所在單元的內(nèi)心為基準(zhǔn)點，給出如下所示的篩選條件的關(guān)系式：

(1)

式中，l為當(dāng)前四面體單元的內(nèi)心與受力點所在單元的內(nèi)心之間的距離，a為篩選系數(shù)，R為受力單元的外接圓半徑。

判斷當(dāng)前的四面體單元是否滿足篩選條件關(guān)系式（1），如果是，將當(dāng)前四面體單元列入?yún)⑴c計算單元的集合，判斷首先從受力點所在四面體單元的相鄰四面體單元開始進(jìn)行，然后再對每個相鄰四面體單元相鄰的四面體單元進(jìn)行判斷。如圖2所示，受力點所在的四面體單元包含i、j、k和p四個點，并與相鄰的四面體單元共用由i、k和p組成的三角形單元，點n1到n2的距離即為公式(1)中的，并且滿足時，可判斷該相鄰四面體單元為參與計算的單元。參與計算單元集合的判斷以受力點所在四面體單元為中心，由里到外進(jìn)行，直到不再滿足篩選條件關(guān)系式（1）為止。篩選參與計算的單元時，根據(jù)篩選判斷的先后順序，將參與計算的單元集合的四面體單元序號和節(jié)點序號進(jìn)行重新排序。當(dāng)受力點所在四面體單元個數(shù)大于1時，依次進(jìn)行受力點所在四面體單元個數(shù)次的循環(huán)進(jìn)行判斷，并統(tǒng)一封裝到參與計算單元的集合。每個四面體單元的相鄰四面體單元，通過步驟1）中生成的數(shù)組A和數(shù)組B進(jìn)行查找，即通過數(shù)組A尋找當(dāng)前四面單元所包含的四個三角形面單元，再根據(jù)得出的四個三角形面單元尋找與當(dāng)前面四面體單元相鄰的四個四面體單元，從而實現(xiàn)相鄰單元的快速查找，并依次進(jìn)行篩選條件的判斷。

步驟3）：創(chuàng)建數(shù)組C，將步驟2）中生成的參與計算單元集合的新四面體單元序號和原模型中的四面體單元序號相連，創(chuàng)建數(shù)組D，將步驟2）中生成的參與計算單元集合的新節(jié)點序號和原模型中的節(jié)點序號相連，能夠?qū)崿F(xiàn)快速查詢篩選出的參與計算的四面體單元新序號和節(jié)點新序號與原模型對應(yīng)的原序號，并能夠保證實時更新發(fā)生形變的節(jié)點坐標(biāo)。

步驟4）：基于步驟2）所生成的參與計算的單元集合，并結(jié)合步驟3）中所生成的數(shù)組C和數(shù)組D，封裝整體剛度矩陣，實現(xiàn)根據(jù)施加的外加載荷和邊界條件，只對軟組織模型中所篩選出的參與計算單元的集合進(jìn)行形變計算。

具體實施方式二：本實施方式所述的一種軟組織有限元模型的加速計算方法，所述的步驟1）中的數(shù)組A為m行4列的二維數(shù)組，數(shù)據(jù)類型為int型，m為軟組織模型所包含的四面網(wǎng)格單元的個數(shù)，每個行下標(biāo)代表軟組織模型的四面體網(wǎng)格單元的序號，每行的四個列下標(biāo)所對應(yīng)的元素分別存儲行下標(biāo)對應(yīng)四面體單元所包含的四個三角形單元的序號。例如，A[2][2]存儲序號為2（假設(shè)四面體單元序號從0開始）的四面體單元中第三個三角形單元的序號。數(shù)組B為n行2列的二維數(shù)組，數(shù)據(jù)類型為int型，n為軟組織模型所包含的三角形單元的個數(shù)，每個行下標(biāo)代表軟組織模型的三角形網(wǎng)格單元的序號，每行的兩個列下標(biāo)所對應(yīng)的元素分別存儲共用當(dāng)前三角形單元的兩個四面體單元的序號，如果當(dāng)前三角形單元為軟組織模型邊界面上的單元，則第一個元素存儲包含當(dāng)前三角形單元的四面體單元的序號，并且第二個元素為0。例如，B[2][1]存儲共用序號為2（假設(shè)三角形單元序號從0開始）的三角形單元的兩個相鄰四面體單元中，第2個四面體單元的序號，如果該三角形單元為軟組織模型的邊界單元，則B[2][1]存儲的值為-1。

具體實施方式三：本實施方式所述的一種軟組織有限元模型的加速計算方法，所述步驟3)中的數(shù)組C為int型的一維數(shù)組，數(shù)組C的下標(biāo)表示原模型的四面體單元序號，每個參與計算單元的新四面體單元序號分別存儲到對應(yīng)的原模型四面體單元序號下標(biāo)的數(shù)組C的元素中，而原模型四面體單元序號中，非參與計算的單元下標(biāo)的對應(yīng)元素值賦為-1。例如，參與計算單元集合中序號為2的四面體單元在原模型中的序號為7（假設(shè)四面體單元序號從0開始），則C[7]賦值為2，如果原模型中序號為7的四面體單元不是參與計算單元集合中的四面體單元，則C[7]賦值為0。數(shù)組D為int型的一維數(shù)組，數(shù)組D的下標(biāo)表示原模型的節(jié)點序號，每個參與計算單元的新節(jié)點序號分別存儲到對應(yīng)的原模型節(jié)點序號下標(biāo)的數(shù)組D的元素中，而原模型節(jié)點序號中，非參與計算節(jié)點下標(biāo)的元素值賦為0。例如，參與計算單元集合中序號為2的節(jié)點在原模型中的序號為7（假設(shè)節(jié)點序號從0開始），則D[7]賦值為2，如果原模型中序號為7的節(jié)點不在參與計算單元集合中的四面體單元上，則D[7]賦值為-1。

具體實施方式四：本實施方式所述的一種軟組織有限元模型的加速計算方法，所述步驟4）中，如圖3所示，只對軟組織的有限元模型中參與計算單元集合，即圖3中的灰色區(qū)域進(jìn)行計算。軟組織的模型采用指數(shù)多項式混合形式的超彈性模型，模型的表達(dá)式如下：

(2)

式中，W為應(yīng)變能密度函數(shù)，C₁、β為模型的材料常數(shù)，I₁和I₂為應(yīng)變張量的主不變量；

基于軟組織的指數(shù)多項式混合形式的超彈性模型建立平衡方程，然后施加邊界條件和外加載荷，可進(jìn)行形變計算。在進(jìn)行形變計算的過程中，平衡方程采用完全拉格朗日法，進(jìn)行分步載荷求解，每一載荷步的非線性方程組采用修正的牛頓迭代法進(jìn)行求解，而在求解非線性方程組的過程中產(chǎn)生的線性方程組可采用CG法進(jìn)行求解。在軟組織有限元計算的過程中，整體剛度矩陣采用GPU計算，非線性方程組采用CPU進(jìn)行求解，來進(jìn)行進(jìn)一步地提高計算效率。