石躍祥 胡 維

1(湘潭大學(xué)信息工程學(xué)院 湖南 湘潭 411105)2(LED照明驅(qū)動(dòng)與控制應(yīng)用工程技術(shù)研究中心 貴州 銅仁 554300)

0 引 言

基于有限元方法的正交各向異性形變體的仿真廣泛應(yīng)用在影視動(dòng)畫(huà)、游戲特效、虛擬現(xiàn)實(shí)等領(lǐng)域中，這些變形效果很大程度上依賴(lài)于材料的本構(gòu)模型[1]，即可以用來(lái)表述形變體材料中的應(yīng)力與應(yīng)變之間的函數(shù)關(guān)系。文獻(xiàn)[2]列出了幾種比較常見(jiàn)的本構(gòu)模型，比如：線(xiàn)性模型、共旋線(xiàn)性模型以及非線(xiàn)性模型等。有了這些材料模型之后，我們可以進(jìn)一步進(jìn)行深入研究，嘗試調(diào)節(jié)材料的一些相關(guān)參數(shù)來(lái)設(shè)計(jì)出一些不同的形變體。但是僅通過(guò)調(diào)這些參數(shù)來(lái)得到更加真實(shí)的形變效果是很困難的。事實(shí)上，這些特殊的標(biāo)準(zhǔn)材料是整個(gè)各向同性材料空間的一小部分，而且整個(gè)各向異性材料空間更加廣泛。

為了獲得更加豐富有趣的形變體形變效果，一些方法已經(jīng)被提出。作為連續(xù)介質(zhì)力學(xué)的基本方程，本構(gòu)方程表示了一種給定材料的物理特性，這其中包括它引起的與之相關(guān)的材料變形反應(yīng)(例如力、應(yīng)力、能量)[3-5]。然而，這些方法由于Valanis-Landel應(yīng)變能密度模型的限制，都是集中在類(lèi)橡膠材料或者肌肉材料之上。在之后的研究工作中，文獻(xiàn)[6]改進(jìn)了Valanis-Landel應(yīng)變能密度函數(shù)模型。基于此模型，本文提出了穩(wěn)定的帶約束的超彈性材料設(shè)計(jì)方法。但是與文獻(xiàn)[6]不同的是：本文更關(guān)注穩(wěn)定的材料編輯。

材料編輯的一個(gè)主要困難是如何編輯穩(wěn)定的材料，因?yàn)闊o(wú)效的材料會(huì)導(dǎo)致仿真不穩(wěn)定。因此，本文給出設(shè)計(jì)正交各向異性材料的穩(wěn)定性條件，并詳細(xì)討論條件。設(shè)計(jì)穩(wěn)定的各向異性材料比較困難，故提出了一種創(chuàng)新的各向異性材料能量，基于這種能量模型，可以很容易地設(shè)計(jì)具有穩(wěn)定性和明顯各向異性的各向異性材料。

1 相關(guān)工作

形變體仿真在計(jì)算機(jī)圖形學(xué)領(lǐng)域當(dāng)中一直是一個(gè)相當(dāng)重要的研究課題。幾種基于物理的方法被用于變形固體的仿真[7]，其中有限元方法被廣泛使用。文獻(xiàn)[2]詳細(xì)介紹了經(jīng)典的有限元方法，其中最基本的有限元方法模型是線(xiàn)彈性有限元模型。該模型很容易設(shè)計(jì)出來(lái)，因此廣泛應(yīng)用于計(jì)算機(jī)圖形學(xué)領(lǐng)域的初期仿真。考慮到線(xiàn)性有限元模型在仿真過(guò)程中是用柯西小應(yīng)變來(lái)度量形變，因此導(dǎo)致該模型僅適用于彈性體的小變形仿真。在很多形變情況下，彈性體較大的形變都是由于旋轉(zhuǎn)導(dǎo)致的，因此文獻(xiàn)[8]提出了共旋線(xiàn)性有限元模型，通過(guò)極分解，將彈性體局部旋轉(zhuǎn)獨(dú)立出來(lái)。然而，當(dāng)彈性體處于較大拉力或者壓力作用從而產(chǎn)生了較大拉伸形變時(shí)，共旋線(xiàn)性有限元仍然不能正確處理大形變，解決該問(wèn)題的一種常見(jiàn)方法是使用二次的格林應(yīng)變，也就是非線(xiàn)性有限元模型。文獻(xiàn)[9]講述了幾種比較常見(jiàn)的非線(xiàn)性有限元模型。文獻(xiàn)[10]在布料動(dòng)畫(huà)仿真中第一次使用非線(xiàn)性有限元方法。對(duì)于形變體，特別是對(duì)于軟物質(zhì)的形變仿真，在形變過(guò)程中單元很可能發(fā)生翻轉(zhuǎn)的情況。然而單元的翻轉(zhuǎn)通常是一種不正確的狀態(tài)，應(yīng)力計(jì)算常常會(huì)導(dǎo)致計(jì)算有誤。文獻(xiàn)[11]針對(duì)該問(wèn)題進(jìn)行了深入研究，最終提出了一種可翻轉(zhuǎn)有限元方法，用以處理那些在仿真時(shí)出現(xiàn)翻轉(zhuǎn)的單元。文獻(xiàn)[12]使用隱式時(shí)間積分格式在可翻轉(zhuǎn)有限元方法中求解形變靜態(tài)平衡位置，同時(shí)為了使仿真穩(wěn)定以及滿(mǎn)足剛度矩陣是正定的，進(jìn)一步修改了牛頓-拉夫森迭代算法?；谖墨I(xiàn)[11-12]，一些學(xué)者提出了多種改進(jìn)效率的方法，比如降維仿真[13]、材料粗糙化[14]、多重網(wǎng)格法[15]等，這些方法幾乎都是在形變不變量空間下進(jìn)行計(jì)算。本文使用的有限元模型是在形變拉伸空間下進(jìn)行仿真。

很多學(xué)者一直以來(lái)專(zhuān)心于研究有關(guān)各向異性材料方面的問(wèn)題，例如：文獻(xiàn)[16]介紹了橫向各向同性超彈性材料,即形變體在垂直軸線(xiàn)的任一方向上所表現(xiàn)出來(lái)的性質(zhì)是相同的。在非線(xiàn)性有限元方法的基礎(chǔ)之上，文獻(xiàn)[17]提出了一種橫向各向同性材料的算法。文獻(xiàn)[18]為了使橫向各向同性材料在仿真過(guò)程中變得更加快速，在GPU上完成了有關(guān)形變體仿真算法的研究。文獻(xiàn)[11]研究了怎樣能夠在橫向各向同性材料中應(yīng)用可翻轉(zhuǎn)有限元模型。之前的仿真大部分集中于各向同性材料，文獻(xiàn)[19-20]借鑒了在工程力學(xué)中用來(lái)編輯正交各向異性材料的相關(guān)思想，深入研究材料的相關(guān)性質(zhì)以及限制一些物理參數(shù)，最終讓仿真出來(lái)的形變體能夠很好的保持穩(wěn)定性，并在此基礎(chǔ)之上，研究出了適用于一般線(xiàn)性各向異性材料的相關(guān)穩(wěn)定性條件。本文與文獻(xiàn)[19-20]的目標(biāo)都是設(shè)計(jì)穩(wěn)定的材料，但不同的是，本文的目的是設(shè)計(jì)非線(xiàn)性材料。

材料設(shè)計(jì)是一個(gè)有趣的熱門(mén)研究課題,多樣的真實(shí)的彈性體材料一直為研究學(xué)者所追求。目前為止已經(jīng)提出了很多種關(guān)于材料設(shè)計(jì)的方法。在工程物理學(xué)科中，基于物理的材料編輯方法(physically-based modeling method)最為常見(jiàn)。該方法是通過(guò)直接設(shè)計(jì)材料的本構(gòu)函數(shù)，并調(diào)整參數(shù)，使得本構(gòu)函數(shù)曲線(xiàn)擬合實(shí)驗(yàn)數(shù)據(jù)。由于文獻(xiàn)[21]提出的應(yīng)變能密度函數(shù)模型是分離形勢(shì)的，而且數(shù)學(xué)推導(dǎo)很?chē)?yán)格，物理解釋也很直觀，因此該模型被研究學(xué)者用來(lái)預(yù)測(cè)形變特征[22-24]。之后文獻(xiàn)[5]在文獻(xiàn)[21]的基礎(chǔ)之上，展開(kāi)了對(duì)不可壓縮的各向同性形變體的研究。而針對(duì)有關(guān)各向異性材料的編輯方法，文獻(xiàn)[3-4]則重點(diǎn)進(jìn)行了一番研究。其中，文獻(xiàn)[3]在其函數(shù)模型中通過(guò)添加各向異性項(xiàng)部分，提出了適合各向異性材料的函數(shù)模型。文獻(xiàn)[25]在計(jì)算機(jī)圖形學(xué)領(lǐng)域中第一次將基于物理的材料設(shè)計(jì)方法引入進(jìn)來(lái)，之后多個(gè)研究工作在其基礎(chǔ)之上進(jìn)行展開(kāi)。如設(shè)計(jì)人體軀干方面的本構(gòu)模型[26]、設(shè)計(jì)手部手指方面的本構(gòu)模型[28]、設(shè)計(jì)人臉?lè)矫娴谋緲?gòu)模型[27]等，但是這些研究進(jìn)展都是針對(duì)某一特定類(lèi)型的材料，并不適用于一般性的材料設(shè)計(jì)。針對(duì)此種情況，文獻(xiàn)[6]運(yùn)用和文獻(xiàn)[5]相似的思路，在設(shè)計(jì)材料的應(yīng)變能密度函數(shù)的時(shí)候參考了文獻(xiàn)[21]的模型來(lái)編輯更加普遍的材料。雖然我們可以基于能量模型和文獻(xiàn)[6]的樣條插值方法上設(shè)計(jì)新的材料，但是設(shè)計(jì)出的材料往往仿真不穩(wěn)定，會(huì)出現(xiàn)固體膨脹、收縮等錯(cuò)誤。仔細(xì)的設(shè)計(jì)可以避免這些問(wèn)題，但是就會(huì)顯得不方便。所以本文詳細(xì)討論了如何更加穩(wěn)定地設(shè)計(jì)材料并且更加方便地展示了穩(wěn)定性條件。除了基于物理的建模方法，我們注意到基于數(shù)據(jù)驅(qū)動(dòng)的材料設(shè)計(jì)方法[29-31]和基于測(cè)量的建模方法[32-34]已經(jīng)在圖形學(xué)方法中大受歡迎，因?yàn)樗鼈兛梢援a(chǎn)生非常逼真的模擬。

2 可分離的應(yīng)變能密度函數(shù)

應(yīng)變能密度函數(shù)用來(lái)表述給定材料的物理特性。通過(guò)設(shè)計(jì)應(yīng)變能密度函數(shù)可以模擬出新的材料。然而，整個(gè)應(yīng)變能密度函數(shù)的空間是很大的，研究這個(gè)空間不是很容易。

2.1 應(yīng)變能密度函數(shù)

基于Valanis-Landel應(yīng)變能模型[21],本文把應(yīng)變能密度函數(shù)ψ作為λ1、λ2和λ3的函數(shù),λi是變形梯度F的奇異值，也表示沿著主軸捕捉變形的主拉伸。此外，對(duì)于各向同性材料，它們?cè)谌魏畏较蛏隙加邢嗤奶匦?。所以?yīng)變能密度函數(shù)ψ(λ1，λ2，λ3)應(yīng)該滿(mǎn)足對(duì)稱(chēng)性的性質(zhì)：

ψ(λ1，λ2，λ3)=ψ(λi1，λi2，λi3)

(1)

式中：i1、i2、i3是(1,2,3)的排列。但在為各向同性材料建模的時(shí)候，很難以滿(mǎn)足對(duì)稱(chēng)性的要求。為了克服ψ的對(duì)稱(chēng)性,文獻(xiàn)[6]擴(kuò)展了Valanis-Landel可分離應(yīng)變能量密度函數(shù),它采取的形式是:

ψ(λ1，λ2，λ3)=f(λ1)+f(λ2)+f(λ3)+g(λ1λ2)+

g(λ1λ3)+g(λ2λ3)+h(λ1λ2λ3)

(2)

式中：f、g、h是一維非線(xiàn)性的標(biāo)量能量密度函數(shù)，分別表示單軸(長(zhǎng)度)、雙軸(面積)和三軸(體積)應(yīng)變。與Valanis-Landel可分離的應(yīng)變能密度函數(shù)不同的是，Valanis-Landel拋棄了g和h,這是因?yàn)槠渲钟谙鹉z材料g=0以及不可壓縮材料h=0。由于f是與材料變形的拉伸有關(guān)，h是與體積有關(guān)，所以本文集中于編輯橡膠材料h=0。當(dāng)設(shè)計(jì)正交各向異性形變體時(shí)，其中的各向同性部分只需要設(shè)計(jì)f和h。另外為了保證ψ的對(duì)稱(chēng)性，這種分離的形式使材料模型編輯更加簡(jiǎn)單和直觀。

2.2 基于拉伸的仿真

之前的有限元模擬方法大多數(shù)集中于不變量空間，即ψ(I1,I2,I3)，比如可翻轉(zhuǎn)有限元法求解大變形[8]，然而我們的能量模型需要仿真拉伸空間ψ(λ1,λ2,λ3)。

(1) 有限元力的計(jì)算。根據(jù)可翻轉(zhuǎn)有限元的工作[8]，首先由式(1)計(jì)算各項(xiàng)同性得第一類(lèi)Piola-Kirchhoff應(yīng)力,然后使用可分離的應(yīng)變能模型計(jì)算：

(3)

(4)

λ2λ3h′(λ1λ2λ3)

(5)

計(jì)算內(nèi)力為:

式中：W是四面體的體積,Dm是參考形狀矩陣。由于W和Dm在仿真中保持不變，可以預(yù)先計(jì)算并存儲(chǔ)它們。

(2) 單元?jiǎng)偠染仃嚨挠?jì)算。整體剛度矩陣是通過(guò)單元?jiǎng)偠染仃嚰色@得的，對(duì)于每一個(gè)形變體中的四面體單元，通常把單元?jiǎng)偠染仃嘖定義為：

(6)

(7)

2.3 超彈性材料的穩(wěn)定性條件

基于分離應(yīng)變能函數(shù)，用戶(hù)能夠方便地設(shè)計(jì)新的材料模型。然而，隨意設(shè)計(jì)函數(shù)可能會(huì)導(dǎo)致模擬不穩(wěn)定。本文首先總結(jié)出了一些關(guān)于設(shè)計(jì)應(yīng)變能密度函數(shù)模型的穩(wěn)定性條件；然后描述了怎樣滿(mǎn)足這些條件來(lái)方便和友好地設(shè)計(jì)材料；最后討論了設(shè)計(jì)正交各向異性材料的方法。

為了能夠設(shè)計(jì)出穩(wěn)定的正交各向異性形變體的應(yīng)變能密度函數(shù)ψ，其應(yīng)該滿(mǎn)足以下約束條件：

(1)ψ是一個(gè)非負(fù)函數(shù)。

(2) 在形變體未形變狀態(tài)下,即λ1=λ2=λ3=1,ψ取得函數(shù)最小值且最小值為0,這個(gè)條件能夠保證形變體在未發(fā)生形變狀態(tài)下的勢(shì)能為0,內(nèi)部無(wú)彈力。

(3)ψ是一個(gè)凹函數(shù)，也就是形變體形變?cè)酱?，其?shì)能越大,內(nèi)部彈力越大，并且這一條件保證了最終計(jì)算的剛度矩陣正定，使得仿真計(jì)算穩(wěn)定[2]。

(4) 在形變體極限形變狀態(tài)下，即λi=0或者λi=∞，ψ的值趨于正無(wú)窮。或者說(shuō)，當(dāng)λi=0時(shí)，形變體內(nèi)部彈力趨于負(fù)無(wú)窮，而當(dāng)λi=∞時(shí)，形變體內(nèi)部彈力趨于正無(wú)窮。即：

(8)

(9)

(10)

(11)

根據(jù)這些穩(wěn)定性條件，我們可以更好地仿真正交各向異性材料。

3 可正交各向異性材料的編輯

各向異性材料在不同的方向上表現(xiàn)出來(lái)的特性是不同的，這與各向同性材料有顯著的區(qū)別。由于各向異性材料具有巨大的材料空間和復(fù)雜性，設(shè)計(jì)出合理的本構(gòu)模型實(shí)屬不易。而正交各向異性材料相比之下就要容易一些，該材料往往會(huì)在三個(gè)正交方向上(稱(chēng)為材料方向m1,m2,m3)展現(xiàn)出不同的剛性特性，并有廣泛的應(yīng)用。文獻(xiàn)[6]提出了有關(guān)正交各向異性材料的分離應(yīng)變能模型?；谶@種能量模型，可以對(duì)正交各向異性材料進(jìn)行建模，但是它很難保持材料穩(wěn)定，即滿(mǎn)足材料穩(wěn)定條件。為了設(shè)計(jì)出更加穩(wěn)定的正交各向異性形變體模型，我們?cè)谖墨I(xiàn)[6]的基礎(chǔ)之上改進(jìn)了正交各向異性材料的可分離應(yīng)變能模型。

ψ=ψiso+ψo(hù)rtho

(12)

(13)

(14)

其中，應(yīng)變能密度函數(shù)ψ包括兩部分：一部分是各向同性能量項(xiàng)，其表現(xiàn)的是各向同性形變體的物理性質(zhì)；另外一部分則是正交各向異性能量項(xiàng)，其反映的是材料的正交性質(zhì)。因?yàn)榉€(wěn)定性條件(1)要求所有的應(yīng)變能密度函數(shù)ψ都為非負(fù)函數(shù)，所以可以簡(jiǎn)單地使式(12)中的所有各向同性能量項(xiàng)ψiso和各向異性能量項(xiàng)ψo(hù)rtho滿(mǎn)足穩(wěn)定條件，從而設(shè)計(jì)出一個(gè)穩(wěn)定的正交各向異性材料。

3.1 各項(xiàng)同性能量項(xiàng)部分的設(shè)計(jì)

3.2 各項(xiàng)異性能量項(xiàng)部分的設(shè)計(jì)

(15)

因?yàn)椋?/p>

(16)

所以：

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

式中：Fi,j表示第j行等于F的第i行、其他位置都為0的矩陣。

4 實(shí) 驗(yàn)

圖1 穩(wěn)定的cube模型的設(shè)計(jì)

(a) 本文設(shè)計(jì)的材料模型(b) stvk材料模型 (c) Neohookean材料模型圖2 形變體turtle對(duì)比實(shí)驗(yàn)

本文通過(guò)對(duì)比實(shí)驗(yàn)來(lái)說(shuō)明材料穩(wěn)定性條件的必要性,以及用本文方法編輯bar模型時(shí)的有效性。不穩(wěn)定的bar模型如圖3所示(973個(gè)頂點(diǎn),2 912個(gè)四面體)。固定bar的頂部除了重力之外沒(méi)有任何外力負(fù)載。在這組對(duì)比示例中,每個(gè)示例的曲線(xiàn)g′和h′都是相同的并滿(mǎn)足條件,材料a滿(mǎn)足所有條件，材料b和c不滿(mǎn)足條件(2)，材料d不滿(mǎn)足條件(3)。調(diào)整f′故意打破能量最小條件(條件(2))和能量凹性條件(條件(3))。應(yīng)變能密度函數(shù)在滿(mǎn)足能量最小條件和能量凹性條件的情況下是被認(rèn)可的。而且,即使bar模型經(jīng)歷大變形，如圖4所示，材料建模方法也很穩(wěn)健。

(a) (b)

(c) (d)圖3 不穩(wěn)定的bar模型(時(shí)間步長(zhǎng)為0.001)

圖4 bar模型的扭曲(時(shí)間步長(zhǎng)為0.001)

(a) 本文正交異性模型(b) 各向同性模型(c) 共旋線(xiàn)性正交異性模型圖5 形變體cube模型對(duì)比實(shí)驗(yàn)

為了說(shuō)明材料的正交各向異性,分別在立方體的m1、m2、m3方向上應(yīng)用初始位移。與各向同性材料和旋轉(zhuǎn)線(xiàn)性正交各向異性材料比較最大拉伸位置時(shí)的狀態(tài),本文的材料確實(shí)表現(xiàn)出非常明顯的正交異性特性,它分別在m1、m2、m3方向顯示出不同的剛度。正交異性dragon(5 803頂點(diǎn),20 140個(gè)四面體)也表現(xiàn)出了良好的穩(wěn)定性和明顯的正交各向異性，如圖6所示。

(a) 各向同性模型(b) 本文正交異性模型(c) 共旋線(xiàn)性正交異性模型圖6 形變體dragon模型對(duì)比實(shí)驗(yàn)

在與各向同性材料的對(duì)比實(shí)驗(yàn)當(dāng)中，設(shè)置了一樣的楊氏模量Eiso和泊松比v。在仿真旋轉(zhuǎn)正交異性恐龍時(shí),E1=Eiso,E2=0.1×Eiso,E3=100×Eiso。與正交各向異性cube模型的例子不同的是,我們分別在x、y和z方向上應(yīng)用初始力而不是初始移位。

此外，本文的非線(xiàn)性正交各向異性和旋轉(zhuǎn)線(xiàn)性正交各向異性立方體幾乎在相同的時(shí)間步時(shí)長(zhǎng)達(dá)到最大拉伸位置。但是,非線(xiàn)性正交各向異性立方體在第57時(shí)間步到達(dá)最壓縮位置,而旋轉(zhuǎn)的線(xiàn)性正交各向異性立方體到達(dá)該位置在第36個(gè)時(shí)間步長(zhǎng)，如圖7所示。

(a) 非線(xiàn)性正交異性模型(b)旋轉(zhuǎn)線(xiàn)性正交異性模型圖7 形變體達(dá)到最壓縮位置時(shí)的對(duì)比實(shí)驗(yàn)

5 結(jié) 語(yǔ)

為了模擬穩(wěn)定材料保證材料的穩(wěn)定性,本文圍繞材料的穩(wěn)定性條件進(jìn)行詳細(xì)討論，并提出了一種更加方便和直觀的材料建模方法。同時(shí)，展示了一個(gè)正交異性材料模型,可用來(lái)設(shè)計(jì)具有穩(wěn)定性以及明顯的正交各向異性特性的正交各向異性材料。本文改進(jìn)了正交各向異性的可分離應(yīng)變能模型材料中的應(yīng)變能密度函數(shù),在此函數(shù)模型中重點(diǎn)研究了正交各向異性能項(xiàng)。為了避免計(jì)算各向異性形變體仿真的復(fù)雜性，本文主要針對(duì)的是正交各向異性形變體的仿真，所以今后將探索更加方便和直觀的方法來(lái)設(shè)計(jì)一般化的各向異性形變體仿真。