虛功原理¶

根據應力與守恆定律中導出的平衡方程式與邊界條件，推導作為連續體力學邊界值問題弱形式的虛功原理。有限元素法的離散化以此弱形式為出發點。本章同時給出當前構形（以 Cauchy 應力與 Almansi 應變線性部分表示）與參考構形（以第二 Piola-Kirchhoff 應力與 Green-Lagrange 應變表示）兩種形式，證明兩者等價後，再確認其歸結至微小變形的形式。

平衡方程式與邊界條件¶

令作用於連續體、每單位質量的體力為 \(\boldsymbol{g}\)，考慮在當前構形中佔據區域 \(\Omega\) 的物體。邊界 \(\Gamma\) 分為位移規定為 \(\bar{\boldsymbol{u}}\) 的幾何邊界 \(\Gamma_B\)，以及表面力規定為 \(\bar{\boldsymbol{t}}\) 的力學邊界 \(\Gamma_t\)，並令 \(\Gamma = \Gamma_B \cup \Gamma_t\)、\(\Gamma_B \cap \Gamma_t = \emptyset\)。對靜態問題，由應力與守恆定律所示的動量守恆式省略慣性項，可得

\[ \nabla_x \cdot \boldsymbol{\sigma} + \rho \boldsymbol{g} = \boldsymbol{0} \quad \text{於} \ \Omega \]

作為平衡方程式。邊界條件為

\[ \boldsymbol{\sigma} \boldsymbol{n} = \bar{\boldsymbol{t}} \quad \text{於} \ \Gamma_t \]

\[ \boldsymbol{u} = \bar{\boldsymbol{u}} \quad \text{於} \ \Gamma_B \]

。以下將虛功原理推導為平衡方程式與力學邊界條件 \(\boldsymbol{\sigma} \boldsymbol{n} = \bar{\boldsymbol{t}}\) 的弱形式。幾何邊界條件 \(\boldsymbol{u} = \bar{\boldsymbol{u}}\) 則藉由測試函數的選擇納入。

當前構形表示的弱形式¶

在弱形式中，分別將未知位移的允許空間與測試函數空間定義為

\[ \mathcal{U} = \{ \boldsymbol{u} \in [H^1(\Omega)]^d \mid \boldsymbol{u} = \bar{\boldsymbol{u}} \ \text{於} \ \Gamma_B \} \]

\[ \mathcal{V} = \{ \delta \boldsymbol{u} \in [H^1(\Omega)]^d \mid \delta \boldsymbol{u} = \boldsymbol{0} \ \text{於} \ \Gamma_B \} \]

其中 \(d\) 為空間維度，\(H^1(\Omega)\) 為函數及其一階弱導數皆平方可積分的 Sobolev 空間，\(\delta\) 表示變分。在當前構形表示中，\(\Omega\) 為變形後的構形；實際數值求解時，會將其拉回參考構形或已知的中間構形處理。

將平衡方程式乘以權重 \(\delta \boldsymbol{u} \in \mathcal{V}\)，並套用 Gauss 散度定理與力學邊界條件，可得當前構形中的虛功原理如下。

\[ \int_{\Omega} \boldsymbol{\sigma} : \delta \boldsymbol{A}_{(L)}\, dv = \int_{\Gamma_t} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma + \int_{\Omega} \delta \boldsymbol{u}^T \rho \boldsymbol{g}\, dv \]

其中 \(\boldsymbol{A}_{(L)}\) 為 Almansi 應變張量的線性部分，定義為

\[ \boldsymbol{A}_{(L)} = \frac{1}{2}\left( \nabla_x \boldsymbol{u} + (\nabla_x \boldsymbol{u})^T \right), \qquad A_{(L)ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \]

其變分為 \(\delta \boldsymbol{A}_{(L)} = \tfrac{1}{2}(\nabla_x \delta \boldsymbol{u} + (\nabla_x \delta \boldsymbol{u})^T)\)。也就是說，要求解 \(\boldsymbol{u} \in \mathcal{U}\)，使其對任意 \(\delta \boldsymbol{u} \in \mathcal{V}\) 均滿足虛功方程式。左邊為內力虛功，右邊為規定表面力與體力所造成的外力虛功。

由於本式寫在變形後（當前構形）的區域上，實際求解時會重新選擇初始構形 \(\Omega_0\)（參考構形）或已知的中間構形作為參考構形，改寫成增量形式後再求解。關於參考構形的具體選擇（Total Lagrange / Updated Lagrange）與增量分解，請參閱增量分析框架。

初始構形表示的弱形式¶

考慮在參考構形中佔據區域 \(\Omega_0\) 的物體，並將其邊界 \(\Gamma_0\) 分為 \(\Gamma_{0B} \cup \Gamma_{0t}\)。將當前構形表示拉回參考構形後，應力－應變共軛對為第二 Piola-Kirchhoff 應力 \(\boldsymbol{S}\) 與 Green-Lagrange 應變 \(\boldsymbol{E}\)。此時，初始構形中的虛功原理為

\[ \int_{\Omega_0} \boldsymbol{S} : \delta \boldsymbol{E}\, dV = \int_{\Gamma_{0t}} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma_0 + \int_{\Omega_0} \delta \boldsymbol{u}^T \rho_0 \boldsymbol{g}\, dV \]

。其中 \(\rho_0\) 為參考構形中的質量密度，依質量守恆式 \(\rho_0 = J\rho\)，與當前構形中的體力表示等價。

當前構形表示與初始構形表示的等價性¶

兩種表示中的內力虛功，透過變形梯度 \(\boldsymbol{F}\) 與體積比 \(J = \det \boldsymbol{F}\) 的轉換彼此一致。亦即，

\[ \int_{\Omega_0} \boldsymbol{S} : \delta \boldsymbol{E}\, dV = \int_{\Omega} \boldsymbol{\sigma} : \delta \boldsymbol{A}_{(L)}\, dv \]

。外力項亦可藉由質量守恆與表面力轉換證明等價。因此，當前構形中的虛功方程式與初始構形中的虛功方程式，是在不同構形下表達同一原理。參照參考構形的求解方法對應 Total Lagrange 法，參照當前構形（前一個已收斂構形）的求解方法則對應 Updated Lagrange 法。

歸結至微小變形¶

在微小變形假設 \(\boldsymbol{F} \approx \boldsymbol{I}\)、\(J \approx 1\) 下，當前構形與參考構形的區別消失，第二 PK 應力與 Cauchy 應力一致（\(\boldsymbol{S} \to \boldsymbol{\sigma}\)），Green-Lagrange 應變與 Almansi 應變線性部分都歸結為微小應變 \(\boldsymbol{\varepsilon}\)。

\[ \boldsymbol{\varepsilon} = \nabla_S \boldsymbol{u} = \frac{1}{2}\left( \nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^T \right), \qquad \varepsilon_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \]

此時，虛功原理歸結為以 Cauchy 應力 \(\boldsymbol{\sigma}\) 與微小應變 \(\boldsymbol{\varepsilon}\) 表示的弱形式

\[ \int_{\Omega} \boldsymbol{\sigma} : \delta \boldsymbol{\varepsilon}\, dV = \int_{\Gamma_t} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma + \int_{\Omega} \delta \boldsymbol{u}^T \rho \boldsymbol{g}\, dV \]

\[ \delta \boldsymbol{u} = \boldsymbol{0} \quad \text{於} \ \Gamma_B \]

。這就是在微小變形、線性彈性靜態分析中直接用於離散化的弱形式（線性彈性靜態分析（入門・附錄）以此形式為出發點，說明從元素剛度 \(\boldsymbol{K}^e\) 的構成到全域方程式 \(\boldsymbol{K}\boldsymbol{U} = \boldsymbol{F}\) 的組裝）。

代入線性彈性本構律 \(\boldsymbol{\sigma} = \boldsymbol{\mathsf{C}} : \boldsymbol{\varepsilon}\)，並以 Voigt 表示法寫成 \(\hat{\sigma} = D\, \hat{\varepsilon}\)，則弱形式變為

\[ \int_{\Omega} \delta \hat{\varepsilon}^T D\, \hat{\varepsilon}\, dV = \int_{\Gamma_t} \delta \boldsymbol{u}^T \bar{\boldsymbol{t}}\, d\Gamma + \int_{\Omega} \delta \boldsymbol{u}^T \rho \boldsymbol{g}\, dV \]

的形式。