高效能元素定式化¶

3 維實體元素的定式化中所示的標準位移型定式化，在套用於近乎不可壓縮的材料或以彎曲為主的薄壁結構時，會出現稱為鎖定（體積鎖定、剪切鎖定）的不合理過度剛性。為避免此問題，FrontISTR 提供 B-bar 法與 F-bar 法（僅替換 B 矩陣或變形梯度的體積部分）、導入內部自由度的非相容元素、將壓力視為獨立未知場的 u-p 混合元素，以及專用於板與梁結構的 MITC 殼元素與梁元素。

本章依元素類型整理這些高效能元素與結構元素的定式化。

B-bar 法¶

將 8 節點線性六面體元素用於近似不可壓縮材料時，單一元素內的應變場可能與體積不變約束相衝突，造成稱為體積鎖定的過度剛性。B-bar 法以元素中心計算的值替換 B 矩陣中對體積膨脹有貢獻的分量，以緩和過度約束 [Hughes1980]。

令在元素中心 \(\boldsymbol{r} = \boldsymbol{0}\) 由形狀函數的空間導數計算所得的 B 矩陣為 \(\bar{\boldsymbol{B}}\)，在積分點 \(\boldsymbol{r}\) 計算的一般 B 矩陣為 \(\boldsymbol{B}(\boldsymbol{r})\)。在節點 \(\alpha\)、自由度 \(i\) 的位移－應變關係中，對體積應變分量 \((\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33})\) 加上

\[ \Delta B_{i\alpha} = \tfrac{1}{3}\bigl(\bar{B}_{i\alpha}(\boldsymbol{0}) - B_{i\alpha}(\boldsymbol{r})\bigr) \]

剪切分量 \((\varepsilon_{12}, \varepsilon_{23}, \varepsilon_{31})\) 則使用一般的 \(\boldsymbol{B}\)。使用由此得到的 B-bar 矩陣組裝元素剛性與內力向量。

FrontISTR 將此定式化專門提供給 8 節點線性六面體（元素 ID 361；參閱元素編號體系），可套用於小變形、Total Lagrange 與 Updated Lagrange。

F-bar 法¶

在有限變形下，體積變化會透過變形梯度 \(\boldsymbol{F}\) 非線性地發揮作用，因此 F-bar 法 [deSouzaNeto1996] 是在變形梯度層級施行與 B-bar 法相當的體積鎖定處理。

令元素中心 \(\boldsymbol{r} = \boldsymbol{0}\) 評估之變形梯度的體積比為 \(J_0 = \det \boldsymbol{F}(\boldsymbol{0})\)，積分點上的變形梯度體積比為 \(J = \det \boldsymbol{F}(\boldsymbol{r})\)，並將積分點上的變形梯度替換為

\[ \bar{\boldsymbol{F}} = \left(\frac{J_0}{J}\right)^{1/3} \boldsymbol{F} \]

如此可得 \(\det \bar{\boldsymbol{F}} = J_0\)，使整個元素內的體積比統一為元素中心值。應力計算與應變－位移矩陣的建立使用替換後的 \(\bar{\boldsymbol{F}}\)；元素切線剛性則包含此替換所產生的附加項。

FrontISTR 將 F-bar 法專門實作於 8 節點線性六面體，可套用於小變形以及 Total Lagrange / Updated Lagrange 的非線性幾何。

非相容元素¶

8 節點線性六面體元素缺少表現彎曲所需的應變模態，因此在以彎曲為主的問題中會出現彎曲鎖定。非相容元素 [Taylor1976] 為補償此缺陷，在元素內部導入額外的位移模態。

除了元素節點位移 \(\boldsymbol{u}^e\) 之外，另在每個元素內導入僅存在於元素內部的非相容模態自由度 \(\boldsymbol{\alpha} \in \mathbb{R}^{9}\)，共 3 個方向 × 3 個模態，並將位移場近似為

\[ \boldsymbol{u}(\boldsymbol{r}) = \sum_{\alpha=1}^{8} N_\alpha^e(\boldsymbol{r})\, \boldsymbol{u}^e_\alpha + \sum_{k=1}^{3} M_k(\boldsymbol{r})\, \boldsymbol{\alpha}_k \]

對自然座標 \(\boldsymbol{r} = (\xi, \eta, \zeta)\)，非相容形狀函數取為 \(M_1 = 1 - \xi^2\)、\(M_2 = 1 - \eta^2\)、\(M_3 = 1 - \zeta^2\)。這些函數不保證元素邊界上的連續性，但會在元素內增加可重現彎曲模態應變的空間。

元素剛性首先組裝成外部自由度－內部自由度的分塊形式

\[ \begin{bmatrix} \boldsymbol{K}_{dd} & \boldsymbol{K}_{d\alpha} \\ \boldsymbol{K}_{\alpha d} & \boldsymbol{K}_{\alpha\alpha} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u}^e \\ d\boldsymbol{\alpha} \end{bmatrix} = \begin{bmatrix} \boldsymbol{F}^e_{\text{ext}} \\ \boldsymbol{0} \end{bmatrix} \]

組裝成此形式後，透過 \(d\boldsymbol{\alpha} = -\boldsymbol{K}_{\alpha\alpha}^{-1}\boldsymbol{K}_{\alpha d}\,d\boldsymbol{u}^e\) 消去內部自由度進行靜態縮約，得到僅含外部自由度的元素剛性

\[ \boldsymbol{K}^e = \boldsymbol{K}_{dd} - \boldsymbol{K}_{d\alpha}\,\boldsymbol{K}_{\alpha\alpha}^{-1}\,\boldsymbol{K}_{\alpha d} \]

並將其交給整體組裝。

FrontISTR 將非相容元素專門實作於 8 節點線性六面體（C3D8IC），可套用於小變形、Total Lagrange 與 Updated Lagrange。

U-P 混合元素¶

B-bar 法與 F-bar 法是在位移型架構內修正體積分量，而 u-p 混合（U-P）元素則採用混合定式化 [Bathe1996]，將壓力 \(\lambda\) 作為獨立於位移的未知場。對近似不可壓縮材料，例如 Poisson 比極接近 0.5 的橡膠狀材料或塑性變形後的金屬，若僅以位移場滿足體積不變約束會造成體積鎖定；將壓力作為獨立變數可緩和此約束。

將應力分離為偏差分量與壓力，

\[ \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\mathrm{dev}} + \lambda\,\boldsymbol{I}, \qquad \boldsymbol{\sigma}_{\mathrm{dev}} = \mathbf{D}_{\mathrm{dev}}\,\boldsymbol{\varepsilon} \]

其中 \(\mathbf{D}_{\mathrm{dev}}\) 為由彈性矩陣中移除與體積模數 \(K\) 成正比之體積部分後得到的偏差彈性矩陣。壓力 \(\lambda\) 與體積應變 \(g = \mathrm{tr}\,\boldsymbol{\varepsilon}\) 透過壓縮率 \(\alpha^{-1} = 1/K\) 的約束條件

\[ g - \alpha^{-1}\lambda = 0 \]

相互關聯。將位移 \(\boldsymbol{u}\) 與壓力 \(\lambda\) 作為未知數離散化後，可得元素聯立方程組

\[ \begin{bmatrix} \mathbf{K}_{uu} & \mathbf{K}_{up} \\ \mathbf{K}_{up}^{T} & \mathbf{K}_{pp} \end{bmatrix} \begin{bmatrix} d\boldsymbol{u} \\ d\lambda \end{bmatrix} = \begin{bmatrix} \boldsymbol{f}_{u} \\ \boldsymbol{f}_{p} \end{bmatrix} \]

其中 \(\mathbf{K}_{uu}\) 包含偏差彈性分量以及（在有限變形時）幾何剛性，\(\mathbf{K}_{up}\) 為耦合體積應變與壓力的矩陣，\(\mathbf{K}_{pp} = -\int \alpha^{-1}\,\boldsymbol{N}_p \boldsymbol{N}_p^{T}\,dV\) 為壓力穩定化項（\(\boldsymbol{N}_p\) 為壓力的形狀函數）。由於壓力自由度封閉於元素內部，

\[ \mathbf{K}_{\mathrm{eff}} = \mathbf{K}_{uu} - \mathbf{K}_{up}\,\mathbf{K}_{pp}^{-1}\,\mathbf{K}_{up}^{T} \]

因此以此形式進行靜態縮約，並將僅含外部（位移）自由度的有效剛性傳給整體組裝。

FrontISTR 將 U-P 元素專門實作於 8 節點線性六面體，每個元素具有 1 個壓力自由度（元素內為常數）。可套用於小變形、Total Lagrange 與 Updated Lagrange。於 Updated Lagrange 中，先使用客觀應力率（Jaumann／Hughes-Winget 型）更新偏差應力，再由靜態縮約求得的值施加壓力 \(\lambda\,\boldsymbol{I}\)。

殼元素¶

薄壁板殼結構使用以 Reissner-Mindlin 板殼理論為基礎的殼元素。低階位移型板殼元素隨厚度變薄，會透過與體積鎖定類似的機制（剪切鎖定）高估橫向剪切應變，使彎曲模態的剛性過度增加。MITC（Mixed Interpolation of Tensorial Components）法 [Dvorkin1984] [Bathe1986] 僅在元素內預先指定的 tying 點重新取樣剪切應變分量，再將取樣值重新插值至元素內，以避免此問題。

MITC 殼元素的節點位於中面上，每個節點具有 6 個自由度：3 個平移分量，加上繞中面法線座標系的 3 個旋轉分量。元素剛性以中面自然座標與厚度方向合計三維的 Gauss 積分計算，板厚 \(h\) 在構成律計算時作為元素屬性給定。

FrontISTR 提供以單一中面層表示的 MITC3（元素 ID 731）、MITC4（741）、MITC9（743），以及沿厚度方向配置兩層節點的層狀殼元素 MITC3-shell361（761，3\(\times\)2 節點、每節點 3 自由度）與 MITC4-shell361（781，4\(\times\)2 節點、每節點 3 自由度）。層狀殼元素的節點自由度僅包含 3 個平移分量，而相當於旋轉自由度的彎曲模態則由兩層配置表現。

梁元素¶

梁架、框架結構等線構件以梁元素離散化。FrontISTR 使用考慮剪切變形的 Timoshenko 梁定式化，將彎曲與剪切都表示為平移與旋轉自由度的函數。

梁元素的每個節點具有 6 個自由度：3 個平移分量，以及繞梁軸與橫向軸的 3 個旋轉分量。元素剛性以沿梁軸的一維數值積分計算。截面積 \(A\) 與彎曲、扭轉方向的截面二次矩 \(I\) 作為梁的截面常數給定為元素屬性，並與材料的 Young 模數 \(E\)、剪切模數 \(G\) 一同構成伸長、彎曲、扭轉與剪切的各剛性係數。

FrontISTR 提供 2 節點直線梁元素（元素 ID 611），以及以 3 個節點表示的 4 節點四面體實體－梁混合元素（641，用於混合自由度）。