基于網(wǎng)格的拉格朗日法可變形模型研究

任川 周國(guó)輝

摘 要：基于物理的可變形模型一直是計(jì)算機(jī)圖形學(xué)與其他相關(guān)學(xué)科綜合研究的重要問題，其主要的模擬技術(shù)有非物理方法和物理方法兩大類。隨著計(jì)算機(jī)硬件的高速發(fā)展和虛擬現(xiàn)實(shí)應(yīng)用需求的不斷提高，基于網(wǎng)格的拉格朗日法作為主要的物理方法之一，日益受到眾多學(xué)者的重視。本文主要從可變形物體的物理建模出發(fā)，首先介紹了可變形物體的幾種常見的建模方法，著重對(duì)基于網(wǎng)格的拉格朗日法中的三種建模方法進(jìn)行了具體的介紹與比較。

關(guān)鍵詞：可變形體；質(zhì)點(diǎn)-彈簧；有限元；邊界元

中圖分類號(hào)：TP391.9 文獻(xiàn)標(biāo)識(shí)號(hào)：A

The Study of Lagrangian Mesh based Methods Deformable Modeling

REN Chuan， ZHOU Guohui

（College of Computer Science and Information Engineering， Harbin Normal University， Harbin 150025，China）

Abstract： Physically based deformable model has been an important issue in the field of Computer Graphics and other related disciplines comprehensive study， the simulation technology is the main non-physical and physical methods in two categories. With the overgrowth of the computer hardware and the improvement of the application， Lagrangian mesh based methods as one of the major physical methods， attracks increasing attention of many scholars. In this paper， starting from the physical modeling of deformable objects， first introduces several common deformable object modeling approach， and focusing on the Lagrangian mesh method， three modeling methods are specific introduction and comparison.

Keywords： Deformable Objects； Mass-Spring； Finite Element； Boundary Element

0 引言

隨著計(jì)算機(jī)圖形學(xué)在眾多學(xué)科綜合研究中的廣泛應(yīng)用，基于物理的可變形模型中形變行為的仿真即已成為時(shí)下研究的熱點(diǎn)之一。該研究主要是根據(jù)物理學(xué)、力學(xué)中的基本規(guī)律和計(jì)算方法，對(duì)物體的變形進(jìn)行模擬，其中涉及到牛頓力學(xué)、彈性力學(xué)，能量守恒原理，差分幾何，數(shù)值計(jì)算等。同時(shí)，該研究還可通過(guò)應(yīng)變能量函數(shù)定量描述物體的質(zhì)量、慣性、重量等特征形變結(jié)果，因而能夠較好地呈現(xiàn)運(yùn)動(dòng)真實(shí)感，更進(jìn)一步地解決了非物理方法中現(xiàn)存的問題。而目前在為大多數(shù)可變形物體進(jìn)行建模時(shí)，物理建模的方式即已成為必然首選。相應(yīng)地，基于物理建模的方法則主要有：歐拉法、基于網(wǎng)格的拉格朗日法、基于無(wú)網(wǎng)格的拉格朗日法、有限狀態(tài)自動(dòng)機(jī)等一系列方法。

1可變形體建模

1.1 可變形體的概念

所謂可變形體是相對(duì)剛體來(lái)說(shuō)的，剛體具有物體內(nèi)部任意兩點(diǎn)間距離保持不變的特點(diǎn)，這就意味著，剛體不會(huì)產(chǎn)生形變；而可變形體，在物體的內(nèi)力作用下就會(huì)產(chǎn)生形變，即內(nèi)部任意兩點(diǎn)間的距離無(wú)法保持不變。

1.2 建模的方法分析

對(duì)虛擬物體實(shí)施研究時(shí)，首先要對(duì)其進(jìn)行虛擬結(jié)構(gòu)表示，即建立模型。物體的建模技術(shù)主要包括：幾何形狀建模、動(dòng)力學(xué)建模、物理構(gòu)造建模和行為建模。而在虛擬環(huán)境中的建模方法則有基于幾何的建模和基于物理的建模以及兩者混合的建模方法。具體地，幾何模型是用來(lái)表示物體的幾何信息的，包含了幾何信息的數(shù)據(jù)結(jié)構(gòu)和相應(yīng)的算法等重要信息。并且，由于其建模技術(shù)操作簡(jiǎn)單，效率更高等優(yōu)點(diǎn)，現(xiàn)已得到了廣泛的應(yīng)用。而基于物理的建模則遵循物理規(guī)律，相對(duì)于簡(jiǎn)單的幾何模型通過(guò)添加物理學(xué)上的約束，相應(yīng)加大了計(jì)算量和數(shù)據(jù)量，從而使模型具有較好的逼真性。研究中進(jìn)一步發(fā)現(xiàn)，在物體之間進(jìn)行交互時(shí)，幾何模型只能通過(guò)已經(jīng)設(shè)定的函數(shù)曲線來(lái)完成模擬交互，遇有其他情景效果表現(xiàn)即不再理想。而物理模型卻因具有物理參數(shù)，當(dāng)情景轉(zhuǎn)換時(shí)，參數(shù)也隨之發(fā)生變化，因此能夠真實(shí)地表達(dá)模擬交互過(guò)程。另外，將幾何模型與物理模型相結(jié)合也是一種優(yōu)秀的建模方法。方法可將物體的大概輪廓用幾何模型定義出來(lái)，再用物理的方法來(lái)精確地定義結(jié)構(gòu)參數(shù)。由此，該方法不僅解決了建模過(guò)程中計(jì)算量大的問題，還可以使物體獲得更好的逼真效果。

2可變形體的物理建模方法

可變形體建模即需要將質(zhì)量、重量、慣性、變形模式等特征與幾何模型和行為法則相結(jié)合，來(lái)完成更真實(shí)的虛擬模型，因此物理建模的方法將可用來(lái)進(jìn)行可變形體仿真。這種方法多是依托實(shí)際問題，建立求解和控制方程，并基于模型的不同特點(diǎn)與要求，選擇建立數(shù)值計(jì)算過(guò)程。迄今為止，基于網(wǎng)格的方法主要有質(zhì)點(diǎn)-彈簧模型、有限元法和邊界元法，以下即是對(duì)這三種方法的詳細(xì)介紹與分析。

2.1 質(zhì)點(diǎn)-彈簧模型

由于質(zhì)點(diǎn)-彈簧系統(tǒng)是模擬所有可變形物體的最簡(jiǎn)單方法，因此廣泛應(yīng)用于基于物理的離散模型。該法在物體上安插很多質(zhì)點(diǎn)，質(zhì)點(diǎn)與質(zhì)點(diǎn)之間具有某種彈性聯(lián)系，通?？捎脧椈蓙?lái)形象地表示這種彈性聯(lián)系，質(zhì)點(diǎn)間依照一定規(guī)則采用彈簧鏈接，物體隨之將可看成是一張由質(zhì)點(diǎn)構(gòu)成的網(wǎng)，并且質(zhì)點(diǎn)的位置表示物體上某一點(diǎn)的空間位置。其中的質(zhì)點(diǎn)沒有大小，但是物體的質(zhì)量集中在質(zhì)點(diǎn)上，而且質(zhì)量是均勻分布的，同事質(zhì)點(diǎn)之間的彈簧卻是沒有質(zhì)量的。當(dāng)有外力作用于質(zhì)點(diǎn)時(shí)，就會(huì)牽動(dòng)其他相鄰質(zhì)點(diǎn)的運(yùn)動(dòng)來(lái)實(shí)現(xiàn)形變效果。

結(jié)構(gòu)性彈簧、扭曲彈簧和拉伸性彈簧是質(zhì)點(diǎn)-彈簧模型中的三種重要類型，具體如圖1所示。

圖1 質(zhì)點(diǎn)彈簧模型和三種類型彈簧

Fig.1 Mass - spring model and three types of spring

由圖1可見，現(xiàn)對(duì)三種彈簧的作用機(jī)理作如下概述：

（1）將橫向和縱向的各質(zhì)點(diǎn)連接起來(lái)的彈簧叫做結(jié)構(gòu)性彈簧（structural spring），結(jié)構(gòu)性彈簧相當(dāng)于物體的框架，主要固定了物體結(jié)構(gòu)。

（2）將連接對(duì)角線的各質(zhì)點(diǎn)的彈簧稱為扭曲彈簧（shear spring），扭曲彈簧可以控制物體的形變，減少物體扭曲。

（3）連接相鄰間隔一個(gè)質(zhì)點(diǎn)的兩個(gè)質(zhì)點(diǎn)的彈簧為拉伸性彈簧（flexion spring），此彈簧起到了保證物體在形變時(shí)邊緣圓滑的作用。質(zhì)點(diǎn)不僅受到了內(nèi)部彈簧的彈性力，而且還會(huì)受到物體外部的施力的影響。主要的外力就是阻尼力，該力可以在物體發(fā)生形變時(shí)，防止質(zhì)點(diǎn)與質(zhì)點(diǎn)之間的不規(guī)律運(yùn)動(dòng)，因此可用彈簧和阻尼來(lái)連接各個(gè)質(zhì)點(diǎn)。

質(zhì)點(diǎn)-彈簧系統(tǒng)在給定時(shí)間 時(shí)被定義為質(zhì)點(diǎn)（ ）的位置定為 和速度為 ，每個(gè)質(zhì)點(diǎn)所受到的內(nèi)力為 則可通過(guò)計(jì)算相鄰連接的彈簧得到。此外，還有所受外力，如重力、空氣阻力、摩擦力等。每個(gè)質(zhì)點(diǎn)的運(yùn)動(dòng)均滿足牛頓第二定律 。并且對(duì)整個(gè)系統(tǒng)而言：

（1）

這里，M表示質(zhì)點(diǎn)的質(zhì)量矩陣。

物體的運(yùn)動(dòng)情況可以制約物體的形變，所以每個(gè)質(zhì)點(diǎn)的動(dòng)力學(xué)方程都應(yīng)滿足牛頓力學(xué)定律：

（2）

（3）

式中， 表示質(zhì)點(diǎn) 的質(zhì)量； 表示質(zhì)點(diǎn) 的加速度； 為質(zhì)點(diǎn) 的速度； 為粘性系數(shù)； 為質(zhì)點(diǎn)所受的外力； 則為質(zhì)點(diǎn) 和質(zhì)點(diǎn) 之間的彈簧對(duì)質(zhì)點(diǎn) 的應(yīng)力，其中 滿足：

（4）

式中， ， 為彈簧兩端點(diǎn)的位置； 為彈簧的彈性系數(shù)；而 則為彈簧的靜止長(zhǎng)度。

在計(jì)算機(jī)圖形學(xué)中，質(zhì)點(diǎn)-彈簧系統(tǒng)已然不再是傳統(tǒng)意義上的結(jié)構(gòu)模式[1]，也就是經(jīng)由最初的應(yīng)用于靜力學(xué)方程來(lái)模擬人面部的表情，發(fā)展至而今的利用動(dòng)力學(xué)方程來(lái)模擬人體皮膚、器官、肌肉等軟組織的變形，由此而成為虛擬醫(yī)學(xué)治療的關(guān)鍵技術(shù)之一。諸如，文獻(xiàn)[2]使用質(zhì)點(diǎn)彈簧模型構(gòu)造肌肉體，模型中用到了拉格朗日動(dòng)力學(xué)方程和龍格-庫(kù)塔法，而且將肌肉的變形效果模擬得尤為逼真。質(zhì)點(diǎn)-彈簧系統(tǒng)也廣泛應(yīng)用于布料的模擬。Breen在文獻(xiàn)[3]中，通過(guò)對(duì)布料物理性質(zhì)的分析，首度提出了布料的Mass-Spring模型。模型中模擬了布料的拉伸壓縮、沿平面內(nèi)外的剪切彎曲，因而可處理各向同性的力學(xué)行為。

2.2 有限元法

所謂有限元法是一種高效能、常用的計(jì)算方法。對(duì)于復(fù)雜的可變形模型來(lái)說(shuō)，具有高復(fù)雜性的有限元法可以精準(zhǔn)計(jì)算并求解微分方程。因此，這種方法在數(shù)學(xué)界[4]、物理界和工程界[5]均已發(fā)揮至關(guān)重要的作用。有限元法的綜合分析就是將模型分解的單元與整體通過(guò)相互轉(zhuǎn)化，而將連續(xù)域中無(wú)限自由度離散化為有限自由度，同時(shí)將微分方程離散化為代數(shù)方程。

有限元法采納了離散和分片插值的思想，來(lái)處理變形問題。具體地，就是首先依據(jù)求解模型的形狀及實(shí)際問題的物理特點(diǎn)，將連續(xù)的幾何模型域離散化為有限個(gè)某種類型單元，并在單元之間設(shè)置有限個(gè)節(jié)點(diǎn)，用以連接全部的單元，由此構(gòu)成一個(gè)由節(jié)點(diǎn)和單元聯(lián)合并存的集合體。而且在每一單元中，根據(jù)節(jié)點(diǎn)位移條件確定位移模式中的系數(shù)，選取場(chǎng)函數(shù)的節(jié)點(diǎn)值為基本未知量，同時(shí)假設(shè)一個(gè)近似插值函數(shù)來(lái)表示單元中場(chǎng)函數(shù)的分布情況，再通過(guò)求解節(jié)點(diǎn)未知量的有限元方程組和設(shè)置的插值函數(shù)得到節(jié)點(diǎn)的數(shù)值和場(chǎng)函數(shù)。若能實(shí)現(xiàn)對(duì)插值函數(shù)的適當(dāng)選取，即可使函數(shù)的內(nèi)外部與界面滿足于相應(yīng)的條件。其后，又將模型中任意有限元節(jié)點(diǎn)自由度對(duì)應(yīng)整體矩陣分量，也就是將單元矩陣和向量疊加起來(lái)，構(gòu)成一整體的矩陣和向量，由此得到了系統(tǒng)運(yùn)動(dòng)力學(xué)方程。最后，則利用節(jié)點(diǎn)、位移為未知量的線性方程組和彈性力學(xué)的相關(guān)知識(shí)，相應(yīng)解得節(jié)點(diǎn)位移和各個(gè)單元的應(yīng)力，并對(duì)運(yùn)動(dòng)控制方程進(jìn)行求解，由此而進(jìn)一步得到了研究計(jì)算結(jié)果。詳細(xì)實(shí)現(xiàn)過(guò)程如圖2所示。

圖2 有限元法的基本步驟

Fig.2 The basic steps of finite element method

對(duì)不規(guī)則網(wǎng)格的偏微分方程求解主要采取的數(shù)值方法則是有限元法。該法即將可變形物體處理為不規(guī)則網(wǎng)格分割的連續(xù)連接體，并對(duì)偏微分方程進(jìn)行求解。彈性材質(zhì)的偏微分方程可表達(dá)為：

（5）

式中， 表示材質(zhì)本身的密度； 表示所受到的外力，如：重力、碰撞力等； 則表示來(lái)物體內(nèi)部形變所產(chǎn)生的力。

有限元建模就是將對(duì)象在載荷、材料、幾何約束、物理約束等方面的特性進(jìn)行一定的數(shù)學(xué)原理分析以及方法整合后的最終結(jié)果展示。首先需對(duì)問題本身的客觀規(guī)律予以精確描述，并建立優(yōu)化合理的有限元方程。但基于龐大的計(jì)算量考慮，還需保證計(jì)算的效率和精度，此時(shí)就要選取相應(yīng)的算法，并借助計(jì)算機(jī)中的實(shí)用程序?qū)?shù)值進(jìn)行相應(yīng)計(jì)算。接下來(lái)的工作就是有限元建模、數(shù)值計(jì)算、以及結(jié)果的比較處理。

有限元建模主要需對(duì)以下幾個(gè)問題進(jìn)行研究，分別是：?jiǎn)栴}的定義與分析，幾何模型的建立，單元選擇與特性的定義，網(wǎng)格劃分，模型檢查和處理，還包括邊界條件的定義，如圖 3 所示。

圖3 有限元建模的過(guò)程和建模的關(guān)鍵技術(shù)

Fig.3 The process and key technology of finite element modeling

通過(guò)有限元的方法對(duì)變形物體進(jìn)行物理仿真時(shí)，可以取得較好的視覺效果。但是由于過(guò)程中的龐大計(jì)算量，對(duì)硬件的要求通常就會(huì)較高。若將其推演開來(lái)可知，一方面在使用有限元方法模擬彈塑性材料的變形過(guò)程中[6]，需對(duì)四面體網(wǎng)格進(jìn)行有限元建模，并建立線性基本函數(shù)。已有文獻(xiàn)[7]中即提出了應(yīng)變張量，有效計(jì)算了網(wǎng)格阻尼力，而探究得到了物體的有限元顯示求解方法。另一方面，對(duì)于某些易碎，易斷裂的彈塑性材質(zhì)模擬中，可使用通過(guò)線性基本函數(shù)相連的四面體網(wǎng)格[8]，以及文獻(xiàn)[9]提出的應(yīng)變張量，進(jìn)而采用了應(yīng)變率張量來(lái)計(jì)算阻尼力，由此即研發(fā)得到了非線性等式的顯式求解。該方法提供了較好的視覺效果，但卻并不支持實(shí)時(shí)運(yùn)行方式。此時(shí)，若為了支持拓?fù)溥壿嬜兓夷芫S持良好的形狀單元，就可在虛擬節(jié)點(diǎn)算法中復(fù)制原始的單元對(duì)象。

2.3 邊界元法

作為一種常用的數(shù)值分析法，邊界元法不同于有限元法中彈性物體的計(jì)算均在物體內(nèi)部體中完成的模式或方式，而是在邊界面上獲得其真正實(shí)現(xiàn)的。邊界元法先對(duì)物體建立邊界微分方程，并與有限元法相似地進(jìn)行單元?jiǎng)澐?，即網(wǎng)格離散化，而在邊界離散化后，則產(chǎn)生邊界元；其后再將邊界元方程離散為代數(shù)方程，對(duì)其進(jìn)行求解；最終解得邊界積分方程。邊界元法的最主要特點(diǎn)就是可以降低維數(shù)，而且還可將運(yùn)動(dòng)等式的積分形式轉(zhuǎn)化為表面積表達(dá)式。該方法采用了降維，就是將二維問題降為一維問題，以及將三維問題降為二維問題來(lái)提高算法整體運(yùn)行效率。只是，邊界元法的不足之處卻在于：若域內(nèi)存在非線性的方程式，就會(huì)導(dǎo)致邊界方程中的未知量不僅包含邊界上的，還包含區(qū)域內(nèi)的，而此時(shí)，采用邊界元法將不再具有優(yōu)勢(shì)。

在工程中，由于有限元法的局限性而無(wú)法解決某些問題，但邊界元法卻有著高效和精準(zhǔn)的優(yōu)點(diǎn)。若在相同的情況下，邊界元法即比有限元法具有更高的精度，若兩種方法精度相同的情況下，邊界元法則有更高的效率。但邊界元法卻只適用于內(nèi)部材質(zhì)相同的物體，在處理拓?fù)渥兓瘯r(shí)，計(jì)算過(guò)程會(huì)比有限元法更加地繁瑣。如文獻(xiàn)[10]，即采用邊界元法處理表面節(jié)點(diǎn)和邊界值，又利用線彈性函數(shù)和邊界元條件下的彈性動(dòng)力模型，從而成功模擬了物體的形變。

2.4 基于網(wǎng)格的拉格朗日法可變形體物理建模方法的比較

上述即是幾種基本的物理建模方法，而且每種建模方法均有著各自的優(yōu)缺點(diǎn)以及相應(yīng)的適用環(huán)境，對(duì)此即可從算法實(shí)現(xiàn)和計(jì)算量的角度來(lái)考慮選取最為適合的建模方法。經(jīng)由對(duì)以上方法的簡(jiǎn)單對(duì)比，即可得出如下結(jié)論：

第一、質(zhì)點(diǎn)-彈簧模型。建模簡(jiǎn)單、直觀，相對(duì)于有限元法計(jì)算量小，易于求解，且模擬速度較快。但由于精度不如其他方法高，約束處理的未臻自然，因而容易造成彈簧的過(guò)度拉伸。

第二、有限元模型。易于處理復(fù)雜邊界以及高梯度問題，有利于求解高階微分方程，可較好地實(shí)現(xiàn)通用化。但卻并不適合處理計(jì)算量大的流體力學(xué)問題以及求解邊界問題。

第三、邊界元法。求解精度高，可降低維數(shù)問題，在處理無(wú)限域問題中具有出眾優(yōu)勢(shì)，但只可求解均質(zhì)問題。而且由于數(shù)學(xué)工具方法的復(fù)雜度大，而不易達(dá)到通用化。

3.結(jié)束語(yǔ)

基于網(wǎng)格的拉格朗日法是物理建模中的常見方法之一，已經(jīng)涌現(xiàn)了大量卓有成效的研究成果，而且受到專家學(xué)者的廣泛關(guān)注。其建模方法除了上面所述的諸多方法，還包括歐拉法中的基于體的建模和基于面的建模、基于無(wú)網(wǎng)格的拉格朗法中的平滑粒子流體力學(xué)和無(wú)網(wǎng)格法以及有線狀態(tài)自動(dòng)機(jī)中的細(xì)胞自動(dòng)機(jī)等在內(nèi)。在如許大量的建模方法中，卻沒有任何一種模型能夠滿足所有的應(yīng)用需求，因此必須要考慮到模型表示、物理參數(shù)范圍以及實(shí)時(shí)交互式模擬等各種參數(shù)和環(huán)境條件?；诖?，研究人員正將各種方法結(jié)合起來(lái)，運(yùn)用各自的獨(dú)屬優(yōu)勢(shì)來(lái)更好地對(duì)模型進(jìn)行約束，進(jìn)而達(dá)到加強(qiáng)交互控制的顯著效果。

參考文獻(xiàn)：

[1] NEALEN A， MULLER M， KEISER R， et al. Physically based deformable models in computer graphics[J]. Computer Graphics Forum，2006， 25（4）：809-836.

[2] NEDEL L， THALMANN D. Real-time muscle deformations using mass-spring systems[C]//CG Int， 1998： 156-165.

[3] David E. Breen ， Donald H. House ， Michael J. Wozny. Predicting the drape of woven cloth using interacting particles， 1994：365–372.

[4] COURANT R. Variational methods for the solution of problems of equilibrium and vibration. Bull Am Math Soc， 1943，49：l-23.

[5] TRUNER M J， CLOUGH R W， MARTIN H C. Stiffness and deflection analysis of complex structures[J]. J Aeronaut sci， 1956， 23（9）： 505-523.

[6] OBRIEN J F， BARGTEIL A W， HODGINS J K. Graphical modeling and an -mation of ductile fracture[J]. ACM Transactions on Graphics， 2002， 21（3）： 291-294.

[7] MOLINON， BAO Z， FEDKIW R. A virtual node algorithm for changing mesh topology during simulation [J]. ACM Transactions on Graphics， 2004， 23（3）： 385-392.

[8] MULLERM， CHARYPAR D， GROSS M. Particle-based fluid simulation for interactive applications [C]// New York： ACM， 2003： 154-159.

[9] KOESER K， BARTCZAK B， KOCH R. RobustGPU-assisted camera tracking using free-form surface models [J]. Journal of Real Time Image Processing， 2007， 2（2）： 133-147.

[10] JAMES D L， PAID K. Real time simulation of multizone elastokinematic models[C]// IEEE International Conference on Robotics and Automation Washingtom， 2002， 1： 927-932.