基于有限元方法的正交各向異性形變體仿真設計

石躍祥 胡 維

1(湘潭大學信息工程學院 湖南 湘潭 411105)2(LED照明驅動與控制應用工程技術研究中心 貴州 銅仁 554300)

0 引 言

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

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

材料編輯的一個主要困難是如何編輯穩定的材料，因為無效的材料會導致仿真不穩定。因此，本文給出設計正交各向異性材料的穩定性條件，并詳細討論條件。設計穩定的各向異性材料比較困難，故提出了一種創新的各向異性材料能量，基于這種能量模型，可以很容易地設計具有穩定性和明顯各向異性的各向異性材料。

1 相關工作

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

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

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

2 可分離的應變能密度函數

應變能密度函數用來表述給定材料的物理特性。通過設計應變能密度函數可以模擬出新的材料。然而，整個應變能密度函數的空間是很大的，研究這個空間不是很容易。

2.1 應變能密度函數

基于Valanis-Landel應變能模型[21],本文把應變能密度函數ψ作為λ1、λ2和λ3的函數,λi是變形梯度F的奇異值，也表示沿著主軸捕捉變形的主拉伸。此外，對于各向同性材料，它們在任何方向上都有相同的特性。所以應變能密度函數ψ(λ1，λ2，λ3)應該滿足對稱性的性質：

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

(1)

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

ψ(λ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是一維非線性的標量能量密度函數，分別表示單軸(長度)、雙軸(面積)和三軸(體積)應變。與Valanis-Landel可分離的應變能密度函數不同的是，Valanis-Landel拋棄了g和h,這是因為其著手于橡膠材料g=0以及不可壓縮材料h=0。由于f是與材料變形的拉伸有關，h是與體積有關，所以本文集中于編輯橡膠材料h=0。當設計正交各向異性形變體時，其中的各向同性部分只需要設計f和h。另外為了保證ψ的對稱性，這種分離的形式使材料模型編輯更加簡單和直觀。

2.2 基于拉伸的仿真

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

(1) 有限元力的計算。根據可翻轉有限元的工作[8]，首先由式(1)計算各項同性得第一類Piola-Kirchhoff應力,然后使用可分離的應變能模型計算：

(3)

(4)

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

(5)

計算內力為:

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

(2) 單元剛度矩陣的計算。整體剛度矩陣是通過單元剛度矩陣集成獲得的，對于每一個形變體中的四面體單元，通常把單元剛度矩陣K定義為：

(6)

(7)

2.3 超彈性材料的穩定性條件

基于分離應變能函數，用戶能夠方便地設計新的材料模型。然而，隨意設計函數可能會導致模擬不穩定。本文首先總結出了一些關于設計應變能密度函數模型的穩定性條件；然后描述了怎樣滿足這些條件來方便和友好地設計材料；最后討論了設計正交各向異性材料的方法。

為了能夠設計出穩定的正交各向異性形變體的應變能密度函數ψ，其應該滿足以下約束條件：

(1)ψ是一個非負函數。

(2) 在形變體未形變狀態下,即λ1=λ2=λ3=1,ψ取得函數最小值且最小值為0,這個條件能夠保證形變體在未發生形變狀態下的勢能為0,內部無彈力。

(3)ψ是一個凹函數，也就是形變體形變越大，其勢能越大,內部彈力越大，并且這一條件保證了最終計算的剛度矩陣正定，使得仿真計算穩定[2]。

(4) 在形變體極限形變狀態下，即λi=0或者λi=∞，ψ的值趨于正無窮。或者說，當λi=0時，形變體內部彈力趨于負無窮，而當λi=∞時，形變體內部彈力趨于正無窮。即：

(8)

(9)

(10)

(11)

根據這些穩定性條件，我們可以更好地仿真正交各向異性材料。

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

各向異性材料在不同的方向上表現出來的特性是不同的，這與各向同性材料有顯著的區別。由于各向異性材料具有巨大的材料空間和復雜性，設計出合理的本構模型實屬不易。而正交各向異性材料相比之下就要容易一些，該材料往往會在三個正交方向上(稱為材料方向m1,m2,m3)展現出不同的剛性特性，并有廣泛的應用。文獻[6]提出了有關正交各向異性材料的分離應變能模型。基于這種能量模型，可以對正交各向異性材料進行建模，但是它很難保持材料穩定，即滿足材料穩定條件。為了設計出更加穩定的正交各向異性形變體模型，我們在文獻[6]的基礎之上改進了正交各向異性材料的可分離應變能模型。

ψ=ψiso+ψortho

(12)

(13)

(14)

其中，應變能密度函數ψ包括兩部分：一部分是各向同性能量項，其表現的是各向同性形變體的物理性質；另外一部分則是正交各向異性能量項，其反映的是材料的正交性質。因為穩定性條件(1)要求所有的應變能密度函數ψ都為非負函數，所以可以簡單地使式(12)中的所有各向同性能量項ψiso和各向異性能量項ψortho滿足穩定條件，從而設計出一個穩定的正交各向異性材料。

3.1 各項同性能量項部分的設計

3.2 各項異性能量項部分的設計

(15)

因為：

(16)

所以：

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

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

4 實 驗

圖1 穩定的cube模型的設計

(a) 本文設計的材料模型(b) stvk材料模型 (c) Neohookean材料模型圖2 形變體turtle對比實驗

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

(a) (b)

(c) (d)圖3 不穩定的bar模型(時間步長為0.001)

圖4 bar模型的扭曲(時間步長為0.001)

(a) 本文正交異性模型(b) 各向同性模型(c) 共旋線性正交異性模型圖5 形變體cube模型對比實驗

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

(a) 各向同性模型(b) 本文正交異性模型(c) 共旋線性正交異性模型圖6 形變體dragon模型對比實驗

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

此外，本文的非線性正交各向異性和旋轉線性正交各向異性立方體幾乎在相同的時間步時長達到最大拉伸位置。但是,非線性正交各向異性立方體在第57時間步到達最壓縮位置,而旋轉的線性正交各向異性立方體到達該位置在第36個時間步長，如圖7所示。

(a) 非線性正交異性模型(b)旋轉線性正交異性模型圖7 形變體達到最壓縮位置時的對比實驗

5 結 語

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