DynEarthSol2D 計算步驟

DynEarthSol2D 數值計算可計算連續體中每一個節點在不同時階(Time Step)之 控制方程式(Governing Equations)，時階為程式反覆運算之次數，依 Choi et al. (2013) 將軟體數值計算可分為六個步驟，以下將依序介紹之，方程式中粗體字代表為向量 或張量，需注意方向正負值。

一、 運動方程式(Equation of Motion)

DynEarthSol2D 模擬之運動方式遵循線性動量平衡，如下列方程式：

𝜌𝒖̇ = 𝛻 ∙ 𝝈 + 𝜌𝒈 (式 3.20) 式3.20 中，𝜌為物質密度；𝒖為速度向量；𝝈為應力張量；𝒈為重力加速度；𝛻表示 為空間梯度；𝛻 ∙可表示為散度運算子(divergence operator)，物體運動依據拉格朗法 (Lagrangian formulation)，以三角形非結構性網格(unstructured mesh)利用上式離散 化，並考慮邊界條件影響，而為了求得節點a 的總力(total force, 𝒇_𝑎)，需將式 3.20 乘上權重函數並做積分，可得到各節點a 的加速度方程式 aa如下式：

𝑚_𝑎𝒂_𝑎 = 𝒇_𝑎 = 𝒇_𝑎^𝑖𝑛𝑡+ 𝒇_𝑎^𝑏𝑐+ 𝒇_𝑎^𝑒𝑥𝑡 (式 3.21) 其中節點的質量𝑚_𝑎為：

𝑚_𝑎 = ∑ (∫ 𝑁_𝑒^𝑎𝜌_𝑓 邊界條件則依紐曼邊界條件(Neumann boundary condition)，邊界力𝒇_𝑎^𝑏𝑐方程式 如下： 模擬計算，DynEarthSol2D 以動力鬆弛法(dynamic relaxation)解算非線性問題之計 算，可由動態的動量方程式獲得靜態平衡。再利用質量縮放技術(mass scaling)以一

及速度，而由這些數據可計算出每個節點的位置及瞬時速度。以下將詳細介紹動力

此外，利用 CFL 條件(Courant–Friedrichs–Lewy condition)限制時間步階(time step)的大小，對於長期構造模擬來說，需要極多的時間步階才足以達到目標的變形

二、 節點混和離散法(Nodal Mixed Discretization)

在DynEerthSol2D 中，三角網格在不可壓縮的變形時遵循體積恆定(volumetric locking)，但當長期構造演變須考慮到塑性或黏性流動的變形時，則必須利用節點 混和離散法(NMD)來達到反體積恆定(antivolumetric locking)的變形模式。元素 e 的 應變速率𝛆̇_𝑒可由速度計算之，如下式：

𝝐̇_{𝑒,𝑖𝑗}^𝑡+∆𝑡 = 1 質張量(identity tensor)。

節點混和離散法的基本概念為將一群鄰近元素之體積應變速率取平均值，並

將其視為由黏彈 (viscoelastic) 行為及彈塑(elasroplactic) 行為兩個次模型所組合 之外，則須利用Simo and Hughes (2004)所提出的 return-mapping alogorithm 將彈性 應力投影至屈服面上。

而對於Mohr-Coulomb 材質，通常用主應力來描述其剪力破壞的屈服函數，如 下式：

𝑓_𝑠(𝜎₁, 𝜎₃) = 𝜎₁− 𝑁_𝜙𝜎₃ + 2𝐶√𝑁𝜙 (式 3.39)

式 3.39 中，𝜎₁與𝜎₃分別表示最大及最小壓縮主應力，其中伸張為正、壓縮為負， 式3.42 中，𝜓為膨脹角(dilation angle)；而發生張力破裂的可能性則為：

𝑔_𝑡(𝜎₃) = 𝜎₃ − 𝜎_𝑡 (式 3.43) 由於需考慮物質的可塑性 (plasticity)，可將其應變增量∆𝛜定義為：

∆𝛜 = ∆𝝐_𝑒𝑙+ ∆𝝐_𝑝𝑙 (式 3.44) 式3.44 中，∆𝝐_𝑒𝑙及∆𝝐_𝑝𝑙分別代表彈性及塑性應變增量，塑性應變增量∆𝜀_𝑝𝑙垂直於流 勢面 (flow potential surface)，可被表示為：

∆𝝐_𝑝𝑙= 𝛽𝜕𝑔

𝜕𝝈 (式 3.45)

式3.45 中，𝛽為塑性流規模，與作用於屈服面的應力狀態有關。而屈服函數可定義

主應力軸的應力及應變狀態可以𝜎_𝐴 = 𝐸_𝐴𝐵𝜖_𝐴此關係式表示，式中，𝜎_𝐴及𝜖_𝐴分別 代表主應力及主應變；而𝐸代表相對應的彈性模數矩陣 (elastic moduli matrix)，其 定義如下：

形等這種與速率不相依的組構模式(constitutive model)時，因其應變與坐標系有關，

例如物體旋轉時，會因座標改變而造成應變變化，但其本身的應力未改變，因此在 計算應力時，須利用每個階段的應力速率去做積分得到應力，而非直接利用應變計 算應力。DynEarthSol2D 利用 Jaumann stress rate (𝛔̇̌)來表示之：

𝝈̇̌ = 𝛔̇ − 𝝎 ∙ 𝛔 + 𝛔 ∙ 𝛚 (式 3.52)

四、 速度及位移計算(Velocity and Displacement Update)

速度是由阻尼加速度計算之，但仍需考慮邊界條件的速度影響，其定義方程式

五、 熱演化(Modeling Thermal Evolution)

岩石圈的熱演化也是長時間構造模擬的一大重點，其主要採用之熱方程式如下：

ρ𝑐_𝑝𝑇 = 𝑘∇̇ ²𝑇 (式 3.57) 式3.57 中，T 為溫度；𝑐_𝑝及k 分別代表岩石圈的熱容量 (heat capacity) 及熱傳導係 數 (thermal conductivity) ，將上式兩邊皆乘上權重函數並做積分，可得下式：

𝐶_𝑎𝑇̇_𝑎^𝑡+∆𝑡 = − ∑ (𝑘𝐷_𝑎𝑏𝑇_𝑏^𝑡Ω_𝑒) 調整節點，利用edge-flipping 調整網格邊界，此方法能有效運用於 Lagrangian 架構 中較大的變形問題(Braun and Sambridge, 1994)。如圖 3-18，白色實線為原本的邊 界，白色虛線為經edge-flipping 重建後的邊界，粉紅點為新加入的節點。當新的網 格重建完畢後，即可重新計算邊界條件、形狀函數及質量矩陣。

圖3-18 重建網格示意圖(Choi et al., 2013)。 discretization) ，將整個模型的保守場(conservative mapping)簡化成下式：

∫ 𝑥^𝑑Ω

同樣地，下標p 為位於新網格 P 且與 q 無交集的元素，𝑥_𝑝為元素 p 上的變異場 x 的值，而Ω_𝑝∩𝑞為元素 p 與 q 之間重疊的部分(如圖 3-19)。在此利用 Farrel and Maddison(2011)所提出的 supermesh 方法，即利用一個涵蓋元素 q 及其與 p 的交點 的大網格演算出p 以及Ω_𝑝∩𝑞。

圖3-19 重建網格說明圖，實線三角形為新建立之網格，虛線為未經重建之舊網格(Choi et al.,

2013)。