NEWMARK

19

4.1 簡化 St. Venant-Kirchhoff 有限元素法

在諸多有限元素法之中，St. Venant-Kirchhoff (StVK)材質是一個常用在圖學領域的材 質。簡單說明，StVK 材質是指該材質在彈性伸縮是均質性的。而 StVK 材質因為簡單容

20

此為簡化的運動方程式，而 Barbic 和 James[2]所取的有效基底稱為質量正交基底，具有 的特性，其中 是 的單位矩陣。可得簡化的運動方程式:

21

22

23

4.2.1 Courant–Friedrichs–Lewy (CFL) 條件

變分積分器可以有效保持系統能量。但變分積分器還是受限於 Courant–Friedrichs–

24 (explicit-implicit variational integrator 簡稱 IMEXv)[1]，對於低頻透過顯性運算，高頻項 透過隱性運算，來達到有效計算的效果。4.3.1 是保守場 IMEXv 的推導流程，4.3.2 是非

25

顯隱性 Lagrangian 量:

(4.3.2)

此離散 Lagrangian 量目的是把離散位能分成高低頻處理。為了推導出離散變分積分器，

將(4.3.2)帶入(3.3.7)、(3.3.8)可得:

(4.3.3) (4.3.4)

理論上，(4.3.3)、(4.3.4)此二方程式就是 IMEXv。但此二方程式因為 是需 要透過迭代求解，導致此二方程式也需要迭代求解。所以，直接對於此二方程式求解並

26

價於 Stormer-Verlet。所以，此積分器是介於這兩積分器之間的綜合積分器。

然而，Step 1 及 Step 3 是相當類似的。當離散動量不需要紀錄，可以將 Step 1 及 Step 3 合併更新法如下，將(4.3.12)的 取代為 可得:

帶入(4.3.10)可得:

27 守場幾何力學是透過作用量方程式(3.3.3)推導出 Euler-Lagrangian 方程式，得以推導出變 分積分器，而非保守場變分積分器[16]是透過 Lagrange-d' Alembert principle:

(4.3.16)

其所對應的離散 Lagrange-d' Alembert principle 為:

(4.3.17)

3保守場是指在模擬過程不具有外力或是消散力存在

28

其中 是離散 Lagrangian 量、 和 是離散外力，將會在後面做進一步討 論。而類似之前我們推導出離散 Euler-Lagrangian 方程式的方式，我們可以推導出具有 外力的離散 Euler-Lagrangian 方程式:

(4.3.18)

將我們之前所定義的離散顯隱性 Lagrangian 量 取代(4.3.20)、(4.3.21)中的 推導 可得:

29

同樣此二方程式需要透過迭代求解，類次 4.3.1 的方式我們取中繼離散動量:

(4.3.26) (4.3.27)

將(4.3.26)、(4.3.27)代入(4.3.24)、(4.3.25)。而 、 皆由中點法 定義所以兩者是相同的。因此令離散高頻外力為

30

在簡化 St. Venant-Kirchhoff 變形物體

在 4.1 我們已討論簡化 StVK 有限元素法，最後的結論是取得有效基底 ，而基底即廣義

31

(4.4.1)、(4.4.2)需要透過求解非線性方程組。但此方程組的 Jacobian 矩陣於大的時間區 間時容易是一個病態矩陣(ill-condition matrix)，直接透過 Newton-Raphson 求解變得不可 行。所以，我們使用一個類似 Newton-Raphson 的方法作為非線性求解工具。

我們將(4.4.1)中的 由(4.4.2)取代可得:

32

2. 第 12 行，高頻率外力的 Jacobian 矩陣。在我們實作之中，只有施作阻尼力 是依據 Rayleigh 阻尼模型。所對應的阻尼矩陣是

，其中 、 是阻尼係數。

3. 第 7、18 行中的高頻外力 阻尼力。而實際阻尼力經由下列計算:

。其中， 為中點速度可由 、 差分取得，或者是 由離散動量計算， ，因為簡化 StVK 有限元素法，質量矩陣 為 單位矩陣，所以可以取 。

33

34

35

第五章 實驗

本論文所採取的具有施力的 IMEXv(以下簡稱 IMEXv)所具備優點是有效保持系統能 量，且透過將低頻和高頻位能分開運算使得計算速度可以加速，因此我們需要驗證能量 Kernel Library)為線性代數運算函數庫以達到最佳運算效能。在這裡我們考慮兩種積分 器:IMEXv 是根據 4.4.2 所實作，而 NEWMARK 積分器是由 Barbic[24]所提供的原始碼。

除此之外，我們考慮三種基礎設定分別為物體模型、外力條件及環境阻尼:

物體模型 為了驗證我們系統的使用廣度，我們嘗詴對不同可形變模型進行模擬，分別 為:beam、basket、bridge、horse、tower 和 heart，此六種可形變模型，如圖 10 所示。

外力條件 為了測詴系統對不同形式外力，我們測詴三種外力條件分別是衝擊力、固定

36

時間長度施力、週期性施力:

1. 衝擊力，是指給予一個長達 0.003 秒的衝擊力，如圖 10。

2. 固定時間長度施力，是指給予施力為 0.3 秒的方波力，以測詴系統對於非衝擊 力的響應效果，如圖 11。

3. 週期性施力，是指給系統外力大小是隨著 1Hz 的正弦波(sin)變化，以測詴系 統是否會發生共振現象，如圖 12。

環境阻尼 為了測詴系統阻尼環境效應。所以，我們測詴兩個阻尼狀態，一為無阻尼 (Rayleigh 阻尼係數 )；另為具有阻尼(Rayleigh 阻尼係數 )。

37

圖 10 可形變模型既施力示意圖。針對每一個模型，做單點施力。藍色箭頭是施力方向 和大小。

38

圖 11 固定時間長度施力示意圖。給予系統單一大小外力加時間為 0.3 秒。在 0.3 秒之 後即無外力作用，如最右下角所示。時間順序為由左至右，由上至下。

39

圖 12 週期性施力示意圖。藍色箭頭為施力方向和大小會隨著時間而改變，時間順序為 由左至右，由上至下。

40 位時間施力)，因為 IMEXv 和 NEWMARK 的外力 quadrature 係數是不相同的(IMEXv

、NEWMARK )，導致兩者能量似乎是不一樣的。

41

圖 13 衝擊力無環境阻尼能量圖。在此兩者能量是有些許不同是因為兩者的 quadrature 不相同(IMEXv 、NEWMMARK )。

圖 14 固定時間長度施力無環境阻尼能量圖。IMEXv 和 NEWMARK，因為施力長達 100 個 frame，使得 quadurate 係數(IMEXv 、NEWMMARK )的影響下 降，所以能量圖近乎相同。

5.3 實驗二: IMEXv 和 NEWMARK 有環境阻尼下的能量行為比較

實驗二是為了比較 IMEXv 和 NEWMARK 具有阻尼的能量行為。實驗二和實驗一基本 設定都是相同的，不同點是實驗二是具有阻尼力版本。實驗一的最後結論是此二積分器 的能量行為相當類似的，所以此節我們僅舉單一例子討論，其餘能量圖一樣請見附錄。

42

圖 16 是衝擊力具有阻尼力的能量圖，可以發現實際上兩者的能量消散行為是極為類似 的並沒有太大不同。

圖 15 週期性施力無環境阻尼能量圖。此為施力為 1Hz 正弦波力的能量圖。因為給予 外力是一個持續性週期施力，所以系統能量將持續累積。在此施力狀況，IMEXv 和 NEWMARK 的能量圖是幾乎相同的。

圖 16 衝擊力具有環境阻尼能量圖。IMEXv 和 NEWMARK 的能量消散是極類似的。

43

5.4 實驗三:驗證系統能量的保持

本節我們將驗證 IMEXv 是否可以有效保持系統能量，意即系統能量不會隨著長時間模 擬導致發散或者是過度收斂。因此，實驗三將觀察在長時間模擬之下 IMEXv 是否可以 保持能量在一定範圍之內。此實驗使用 horse 為運動模型，給予一個初始的衝擊力如 5.1 所描述。接著在中間過程不給予外力和阻尼力，設定時間區間大小為 0.003 秒，最後我 們讓模擬長達 10000 秒，也就是多達三百萬個 frames。能量圖如圖 17 所示，可以看出 總體能量保持在固定範圍之內，但實際上後 35 秒的平均總體能比前 35 秒增加了 2%；

而高頻能量數值誤差的累積有所衰減，導致後 35 秒比前 35 秒衰退了 28%。

(a) (b)

(c) (d)

圖 17 良好能量行為驗證圖。此為 horse 的能量圖，左邊是前 35 秒、右邊是最後 35 秒，

表示經過了三百萬個 frames。後 35 秒的整體能量比前 35 秒增加了 2%、而高頻能 量消散 28%。其中，施力條件給予衝擊力，中間不具有阻力和外力。時間區間為 0.003 秒模擬時間長達 10000 秒。(a)、(b)圖為系統總能量，(c)、(d)為系統高頻能

44

量。

5.5 實驗四:不同時間區間大小比較

實驗四將驗證時間區間大小是否影響 IMEXv 的能量結果，本實驗使用 beam 作為施作模 型，施力條件為長達 0.3 秒的方波力，在方波力去除之後即無其他外力和阻尼力，在同 條件下除了時間區間大小皆相同，其中設定時間區間大小為 0.01 秒、0.005 秒、0.001 秒。最後，能量圖如圖 18 所示，結果巨觀顯示此三時間區間的能量曲線幾乎是一模一 樣，因此我們推論 IMEXv 巨觀能量行為不會隨著時間區間大小改變。若微觀即如圖 19 所示可發現不同時間區間大小還是有些許不同。

圖 18 不同時間區間大小巨觀能量圖。此為 beam 的能量圖，表示在時間區間大小並不 會影響巨觀能量圖。此圖可以觀察出在我們設定的三個時間區間的能量曲線是非 常類似的。其中我們設定的是時間區間大小為 0.01 秒、0.005 秒、0.001 秒。施力 條件為初始給予長達 0.3 秒的施力，在中間過程不具有其他外力跟阻尼力。

45

圖 19 不同時間區間大小微觀能量圖。此為圖 18 僅觀察 0~1 秒。可以發現三個時間區 間大小雖然相當類似，但還是有些些許不同，其中時間區間越小變化越是快速，

但整體還是相當類似。

5.6 效能比較

在此我們比較 NEWMARK 和 IMEXv 的運算時間，以驗證 IMEXv 是在效能上有所提升。

而 IMEXv 運算速度的提升是藉由將顯性求解低頻運動模式，而非如 NEWMARK 皆採 取隱性求解。所以，IMEXv 的低頻運動模式數量會直接影響到運算速度。如表 1 是時 間區間為 0.003 秒的高頻低頻運動模式數量表。其中，高低頻運動模式數量是根據 4.2.2 的拆解方法所自動決定，而把低頻數量會直接影響到加速倍率如表 2 所示 ，我們可以 發現 IMEXv 確實可以達到加速的目的。其中，以 heart 這個模型達到最大的加速比例是 因為 heart 在 IMEXv 高頻數量只有 3 也就是只有 3 個運動模式透過隱性求解而 NEWMARK 卻需要 30 個運動模式皆透過隱性求解，所以 heart 這個模型 IMEXv 在速度

46

IMEXv 是一個介於 Stormer-Verlet 和隱性中點法的變分積分器，保留了此二變分積分器 的優點，但也意味包含兩個積分器的缺點。

47

48 和 NEWMARK 最大的差別是，IMEXv 會隨著時間區間大小動態區分出顯性和隱性求解 的運動模式並分開處理，而 NEWMARK 僅把所有運動模式透過顯性或是隱性求解。所 以，IMEXv 的收斂速度會因為非線性方程組降維(見 4.4.2)得以變快，最後 IMEXv 可以 達到 1 到 3 倍的加速(見 5.6)相對於 NEWMARK。

49

但實際上，在我們的實作 IMEXv 的本身求解的穩定性並不如 NEWMARK，在我 們實驗一中的 36 組實驗裡面，有三組實驗 IMEXv 並沒有達成完成長達 30 秒的模擬時 間如表 3 所示，其中的原因是非線性求解還是具有本身穩定性的問題，如圖 21 在 6.5 秒之後系統能量變化開始變得不穩定，最後在七秒之後即發散無法有效求解。所以，其 實 IMEXv 在非線性求解是有待加強的部分，可以在未來繼續完成。

6.3 未來展望

我們所提出的低頻高頻拆解方法是透過靜止狀態的自然頻率進行拆解。假如可以提出一 個在執行時間有效率拆解方法。將會使 IMEXv 使用範圍更加廣泛具有更大吸引力。

目前我們將 IMEXv 套用在簡化的 StVK 有限元素法之上。其實 IMEXv 也可以使用 在更多物理模擬系統像是布料、流體、煙模擬之上，但需要一個有效區分出低頻位能和

目前我們將 IMEXv 套用在簡化的 StVK 有限元素法之上。其實 IMEXv 也可以使用 在更多物理模擬系統像是布料、流體、煙模擬之上，但需要一個有效區分出低頻位能和